

*=============================================================================*
* 
*=============================================================================*

	capture log close
	clear
	set more off
	
	set scheme s2color, permanently
	
	global wealthdon 		<SET FOLDER>
	

	global figures 			${wealthdon}/figures
	global temp_estimates 	${wealthdon}/temp_estimates
	global raw 				${wealthdon}/raw
	global projectdatafolder ${wealthdon}/data/projectdatafolder

*=============================================================================*

disp "[$S_DATE $S_TIME] Starting"
*=============================================================================*
*  Open and prepare data for analyses
*=============================================================================*

use ${projectdatafolder}/wdon1_master1.dta
tsset LOPNR year

gen one = 1

destring bokomm, replace

gen double wtax_threshold = .
replace wtax_threshold = 700000		if year == 2010
replace wtax_threshold = 700000		if year == 2011
replace wtax_threshold = 750000		if year == 2012
replace wtax_threshold = 870000		if year == 2013
replace wtax_threshold = 1000000	if year == 2014
replace wtax_threshold = 1200000	if year == 2015
replace wtax_threshold = 1400000	if year == 2016
replace wtax_threshold = 1480000	if year == 2017	
replace wtax_threshold = 1480000	if year == 2018

gen wtax_rate = .
replace wtax_rate = 0.0110 if inrange(year,2010,2013)
replace wtax_rate = 0.0100 if year==2014
replace wtax_rate = 0.0085 if inrange(year,2015,2018)

gen TNW = bel49
gen TVH_pri = prim_lign
gen TVH_sec = sek_lign
gen TVH_cab = fritid
foreach var of varlist TVH_pri TVH_sec TVH_cab {
	replace `var' = 0 if missing(`var') & !missing(ber_brform)
}

gen 	wtax_TVHS_rate = 0.4 if inrange(year,2010,2012)
replace wtax_TVHS_rate = 0.5 if year == 2013
replace wtax_TVHS_rate = 0.6 if year == 2014
replace wtax_TVHS_rate = 0.7 if year == 2015
replace wtax_TVHS_rate = 0.8 if year == 2016
replace wtax_TVHS_rate = 0.9 if inrange(year,2017,2018)
gen wtax_TVHS_rate_2012 = 0.4

gen don_cap = .
replace don_cap = 12000 if inrange(year,2010,2013)
replace don_cap = 16800 if year == 2014
replace don_cap = 20000 if year == 2015
replace don_cap = 25000 if year == 2016
replace don_cap = 30000 if year == 2017
replace don_cap = 40000 if year == 2018
replace don_cap = 50000 if year == 2019

gen 	don_rate = 0.28 if inrange(year,2010,2013)
replace don_rate = 0.27 if inrange(year,2014,2015)
replace don_rate = 0.25 if year==2016
replace don_rate = 0.24 if year==2017
replace don_rate = 0.23 if year==2018

* missing giving data = zeros (nothing reported to tax authorities)
foreach var of varlist gaver_tot gaver_rel gaver_nonrel gaver_medical gaver_social gaver_climate gaver_humrig gaver_domserv gaver_rest {
	replace `var' = 0 if missing(`var') & inrange(year,2012,2018)
}


merge 1:1 LOPNR year using ${raw}/serie_k3371.dta
gen _present2012 = 1 if !missing(LOPNR) & year == 2012
bysort ab (year): egen _tot_present_2012 = total(_present2012)
tab year _tot_pres
drop if _tot_present==0 & _merge==2
drop _merge
gen byte k3371_missing = missing(k3371)
replace k3371 =0 if missing(k3371)

* Household aggregation
bysort hhid year: egen hh_N_wtaxpayers = count(skatt_formue)

sort LOPNR year
by LOPNR (year): egen double hhid_2012 = total(hhid*(year==2012))
by LOPNR (year): egen int hh_N_wtaxpayers_2012 = total(hh_N_wtaxpayers*(year==2012))


gen debt = GJELD
sort hhid_2012 year
foreach var of varlist gaver_aid debt gaver_tot gaver_rel gaver_nonrel gaver_medical gaver_social gaver_climate gaver_humrig gaver_domserv gaver_rest forskning skatt_formue TNW ber_brform TVH_* k3371  {

	by hhid_2012 year: egen hh_`var' = total(`var'),mi
}
ren hh_skatt_formue hh_wtax

gen hh_wtax_dummy = (hh_wtax>0) if !missing(hh_wtax)

gen hh_ETNW = hh_TNW-hh_N_wtaxpayers*wtax_threshold



* Select hh head based on gross wealth and income
bysort hhid year (ber_brform brutto): gen byte hh_head = _n==_N
sort LOPNR year
by LOPNR (year): egen byte hh_head_2012 = total(hh_head * (year==2012) )


gen grossinc = bel_1_1
replace grossinc = brutto if missing(grossinc)

gen 	laborearnings = wyrkinnt  if inrange(year,2015,2018)
replace laborearnings = wyrkinnt  if inlist(year,2010,2011,2012,2014)
replace laborearnings = yrkinnt   if year==2013
replace laborearnings = 0 if missing(laborearnings) & !missing(rentinnt)

gen MNW = ber_brform - debt

gen int age = alder
gen age2 = (alder/50)^2
gen age3 = (alder/50)^3

gen _size = -antpers_i_regstat_famlopenr
bysort LOPNR (_size): gen famsize=_size[1]
sort LOPNR year
drop _size
gen famsize2 = famsize^2

capture drop *gaver_tot_w
gen    gaver_tot_w = min(gaver_tot,   100000) if !missing(gaver_tot)
gen hh_gaver_tot_w = min(hh_gaver_tot,100000) if !missing(hh_gaver_tot)


capture drop *gaver_tot_dummy 
gen gaver_tot_dummy = gaver_tot > 0 if !missing(gaver_tot)
gen hh_gaver_tot_dummy = hh_gaver_tot > 0 if !missing(hh_gaver_tot)

sort hhid_2012 year
by hhid_2012 year: gen hh_grossinc = grossinc[1] + max(grossinc[2],0)
by hhid_2012 year: gen hh_laborearnings = laborearnings[1] + max(laborearnings[2],0)
by hhid_2012 year: egen hh_MNW = total(MNW), mi
by hhid_2012 year: egen hh_rentinnt = total(rentinnt), mi

gen hh_loggrossinc = log(10000*1.02^(year-2012) + hh_grossinc)
gen hh_loglaborearnings = log(10000*1.02^(year-2012) + hh_laborearnings)

sort LOPNR year
by LOPNR (year): egen double hh_TVH_sec_2012 = total(hh_TVH_sec*(year==2012)),mi
by LOPNR (year): egen double hh_TVH_sec_2011 = total(hh_TVH_sec*(year==2011)),mi
by LOPNR (year): egen double hh_TVH_sec_2010 = total(hh_TVH_sec*(year==2010)),mi
by LOPNR (year): egen double hh_TVH_pri_2012 = total(hh_TVH_pri*(year==2012)),mi
by LOPNR (year): egen double hh_TVH_pri_2010 = total(hh_TVH_pri*(year==2010)),mi
by LOPNR (year): egen double hh_TVH_cab_2012 = total(hh_TVH_cab*(year==2012)),mi
by LOPNR (year): egen double hh_TVH_cab_2010 = total(hh_TVH_cab*(year==2010)),mi
foreach var of varlist gaver_aid gaver_tot gaver_rel gaver_nonrel gaver_medical gaver_social gaver_climate gaver_humrig gaver_domserv gaver_rest  {
	by LOPNR (year): egen 	   hh_`var'_2012 = total(hh_`var'*(year==2012)),mi
}

by LOPNR (year): egen 	   hh_forskning_2012 = total(hh_forskning*(year==2012)),mi
by LOPNR (year): egen double hh_TNW_2012 = total(hh_TNW*(year==2012)),mi
by LOPNR (year): egen double hh_TNW_2010 = total(hh_TNW*(year==2010)),mi
by LOPNR (year): egen double hh_debt_2012 = total(hh_debt*(year==2012)),mi
by LOPNR (year): egen double hh_debt_2010 = total(hh_debt*(year==2010)),mi
by LOPNR (year): egen double hh_wtax_2012 = total(hh_wtax*(year==2012)),mi
by LOPNR (year): egen double hh_ETNW_2012 = total(hh_ETNW*(year==2012)),mi
by LOPNR (year): egen double hh_ber_brform_2012 = total(hh_ber_brform*(year==2012)),mi
by LOPNR (year): egen double hh_loggrossinc_2012 = total(hh_loggrossinc*(year==2012)),mi
by LOPNR (year): egen double hh_laborearnings_2012 = total(hh_laborearnings*(year==2012)),mi

gen 	hh_wtax_b12 = min(hh_wtax,1000000*0.01*10) if hh_TNW_2012<=1000000 & !missing(hh_wtax) & !missing(hh_TNW_2012)
replace hh_wtax_b12  = min(hh_wtax,hh_TNW_2012*0.01*10) if hh_TNW_2012>1000000 & !missing(hh_wtax) & !missing(hh_TNW_2012)

gen hh_cTNW_2012 = hh_ber_brform_2012 -(0.75/0.25)*hh_TVH_pri_2012 - (0.6/0.4)*hh_TVH_sec_2012  - hh_debt_2012
replace hh_cTNW_2012 = hh_TNW_2012 if !missing(hh_TNW_2012) & hh_TNW_2012>0


gen hh_cTNW = hh_TNW if !missing(hh_TNW) & hh_TNW>0
replace hh_cTNW = hh_ber_brform_2012 ///
		-(0.75/0.25)*hh_TVH_pri ///
		- (1-wtax_TVHS_rate/wtax_TVHS_rate)*hh_TVH_sec ///
		- hh_debt	if missing(hh_cTNW)
gen hh_cETNW = hh_cTNW-hh_N_wtaxpayers*wtax_threshold

bysort LOPNR (year): egen double age_2012 = total(age*(year==2012)),mi
gen age2_2012 = (age_2012/50)^2
gen age3_2012 = (age_2012/50)^3

by LOPNR (year): egen double famsize_2012 = total(famsize*(year==2012)),mi
gen famsize2_2012 = famsize_2012*famsize_2012

gen hh_logMNW = log(1000*1.02^(year-2012) + hh_MNW) 

*=============================================================================*
* SETTINGS
*=============================================================================*

capture drop _hh_ETNW_2012_bin_100k
gen long _hh_ETNW_2012_bin_100k = 1000 + floor((hh_cTNW_2012-hh_N_wtaxpayers_2012*750000)/100000 )
replace _hh_ETNW_2012_bin_100k  = 1000 + floor((hh_cTNW_2012-hh_N_wtaxpayers_2012*750000)/200000 )*2 if _hh_ETNW_2012_bin_100k >= 1000 + 4
replace _hh_ETNW_2012_bin_100k  = 1000 + floor((hh_cTNW_2012-hh_N_wtaxpayers_2012*750000)/400000 )*4 if _hh_ETNW_2012_bin_100k >= 1000 + 10
replace _hh_ETNW_2012_bin_100k  = 1020 if _hh_ETNW_2012_bin_100k>=1020
replace _hh_ETNW_2012_bin_100k  = 985 if _hh_ETNW_2012_bin_100k<=985
replace _hh_ETNW_2012_bin_100k  = . if missing(hh_cTNW_2012-hh_N_wtaxpayers_2012)

global if_TVHS_TNW = "inrange(hh_TNW_2012,1,2000000) & hh_MVH_prisecM_2012>0 & !missing(hh_MVH_prisecM_2012) "
global if_TVHS_cTNW = "inrange(hh_cTNW_2012/(750000*hh_N_wtaxpayers_2012),-2,4) & hh_MVH_prisecM_2012>0 & !missing(hh_MVH_prisecM_2012) "
global if_TVHS_cTNW_2 = "_hh_ETNW_2012_bin_100k>=995 & hh_cTNW_2012/(750000*hh_N_wtaxpayers_2012) <= 4 & hh_MVH_prisecM_2012>0 & !missing(hh_MVH_prisecM_2012) "
global ifSR = "inrange(gaver_tot-don_cap,-7000,6999)"
global ifSRwide = "inrange(gaver_tot-don_cap,-10000,9999)"
global if_RKD_TNW = "!missing(hh_TVH_pri+hh_TVH_sec+hh_TVH_cab) & TVH_pri+hh_TVH_sec+hh_TVH_cab>0"

*=============================================================================*
* Misc variables
*=============================================================================*

gen G_share_rel = hh_gaver_rel/hh_gaver_tot 

gen FD_log_g = log(1000*1.02^(year-2012) + F.hh_gaver_tot_w) - log(1000*1.02^(year-2012) + hh_gaver_tot_w) 
gen log_g = log(1000*1.02^(year-2012) + hh_gaver_tot_w)
gen D_log_g = D.log_g

gen g_rel = hh_gaver_rel
gen g_nonrel =  hh_gaver_nonrel
gen g_intl = hh_gaver_aid + hh_gaver_climate + hh_gaver_humrig
gen g_dom = hh_gaver_social +hh_gaver_domserv+hh_gaver_medical+ hh_gaver_rest
gen log_g_rel 		= log(1000*1.02^(year-2012) + hh_gaver_rel)
gen log_g_nonrel 	= log(1000*1.02^(year-2012) + hh_gaver_nonrel)
gen log_g_intl 		= log(1000*1.02^(year-2012) + hh_gaver_aid + hh_gaver_climate + hh_gaver_humrig)
gen log_g_dom		= log(1000*1.02^(year-2012) +  hh_gaver_social +hh_gaver_domserv+hh_gaver_medical+ hh_gaver_rest)
gen hh_gaver_intl_2012 = hh_gaver_aid_2012 + hh_gaver_climate_2012 + hh_gaver_humrig_2012

gen hh_wtax_b12_NOK1000 = hh_wtax_b12/1000
gen hh_AWTR = hh_wtax/hh_MNW
replace hh_AWTR = 0 if hh_wtax==0 & !missing(hh_MNW)
replace hh_AWTR = min(max(hh_AWTR,0),wtax_rate) if !missing(hh_AWTR) & !missing(wtax_rate)

gen hh_MWTR = wtax_rate*hh_wtax_dummy
 
gen hh_laborearnings_w = max(min(hh_laborearnings,1500000*1.02*(year-2012)),0) if !missing(hh_laborearnings) /* 1.5M > 99pctile*/

gen hh_k3371_w  = min(max(hh_k3371,0),100000*1.02^(year-2012)) if !missing(hh_k3371) /* 10M max */

gen hh_MVH_secM_2012    = min((hh_TVH_sec_2012/0.4 				)/1000000,10) if !missing(hh_TVH_sec_2012)
gen hh_MVH_secM_2011    = min((hh_TVH_sec_2011/0.4 				)/1000000,10) if !missing(hh_TVH_sec_2011)
gen hh_MVH_secM_2010    = min((hh_TVH_sec_2010/0.4 				)/1000000,10) if !missing(hh_TVH_sec_2010)
gen hh_MVH_cabM_2012 	= min((hh_TVH_cab_2012/0.30 				)/1000000,10) if !missing(hh_TVH_cab_2012)
gen hh_MVH_cabM_2010 	= min((hh_TVH_cab_2010/0.30 				)/1000000,10) if !missing(hh_TVH_cab_2010)
gen hh_MVH_priM_2012 	= min((hh_TVH_pri_2012/0.25 				)/1000000,10) if !missing(hh_TVH_pri_2012)
gen hh_MVH_priM_2010 	= min((hh_TVH_pri_2010/0.25 				)/1000000,10) if !missing(hh_TVH_pri_2010)
gen hh_MVH_prisecM_2012 = hh_MVH_priM_2012 + hh_MVH_secM_2012
gen hh_MVH_cabsecM_2012 = hh_MVH_cabM_2012 + hh_MVH_secM_2012
gen hh_MVH_cabsecM_2010 = hh_MVH_cabM_2010 + hh_MVH_secM_2010
gen hh_MVH_pricabM_2012 = hh_MVH_cabM_2012 + hh_MVH_priM_2012
gen hh_MVH_priseccabM_2012 = hh_MVH_priM_2012 + hh_MVH_cabM_2012 + hh_MVH_secM_2012
gen hh_MVH_priseccabM_2010 = hh_MVH_priM_2010 + hh_MVH_cabM_2010 + hh_MVH_secM_2010


by LOPNR (year): egen int _N_posgiving_pre = total( (hh_k3371>0)*!missing(hh_k3371)*inrange(year,2006,2012) )
gen byte hh_pre_giver = .
replace hh_pre_giver = 1 if _N_posgiving_pre == 7
replace hh_pre_giver = 0 if _N_posgiving_pre == 0

gen hh_MVH_secM_adj_2012= min(hh_MVH_secM_2012*wtax_TVHS_rate,5) if !missing(hh_MVH_secM_2012)

gen hh_MVH_prisecM_2012_sq 	= hh_MVH_prisecM_2012^2
gen hh_MVH_cabsecM_2012_sq 	= hh_MVH_cabsecM_2012^2
gen hh_MVH_secM_2012_sq 	= hh_MVH_secM_2012^2
gen hh_TNW_2012_sq		= (hh_TNW_2012/1000000)^2
gen hh_MVH_secM_adj_2012_sq	= hh_MVH_secM_adj_2012^2

gen hh_TNW_2012_cu		= (hh_TNW_2012/1000000)^3

gen hh_MVH_sec_dummy_2012 = hh_MVH_secM_2012 > 0 if !missing(hh_MVH_secM_2012)
gen hh_MVH_cabsec_dummy_2012 = hh_MVH_cabsecM_2012 > 0 if !missing(hh_MVH_cabsecM_2012)

gen dt = year-2012

gen Post = year>=2013

capture drop not2012
gen not2012 = year!=2012

* Robustness check variables
gen hh_gaver_tot_w_R1 = hh_gaver_tot_w  
gen hh_wtax_b12_R1 = hh_wtax_b12
gen D2012_2010_hh_MVH_secM = hh_MVH_secM_2012 - hh_MVH_secM_2010


gen hh_MVH_prisec = TVH_pri/0.25 + TVH_sec/wtax_TVHS_rate
gen hh_logMVH_prisec = log(hh_MVH_prisec)

*-------------------------------------------------*
* Interaction variables (e.g., TNW_{2012} bins)
*-------------------------------------------------*

capture drop _hh_ETNW_2012_bin_wide
gen long _hh_ETNW_2012_bin_wide = 1000 + floor((hh_cTNW_2012-hh_N_wtaxpayers_2012*750000)/200000 )*2
replace _hh_ETNW_2012_bin_wide  = 1000 + floor((hh_cTNW_2012-hh_N_wtaxpayers_2012*750000)/400000 )*4 if _hh_ETNW_2012_bin_wide >= 1000 + 4
replace _hh_ETNW_2012_bin_wide  = 1000 + floor((hh_cTNW_2012-hh_N_wtaxpayers_2012*750000)/800000 )*8 if _hh_ETNW_2012_bin_wide >= 1000 + 8
replace _hh_ETNW_2012_bin_wide  = 1020 if _hh_ETNW_2012_bin_wide>=1020
replace _hh_ETNW_2012_bin_wide  = 985 if _hh_ETNW_2012_bin_wide<=985
replace _hh_ETNW_2012_bin_wide  = . if missing(hh_cTNW_2012-hh_N_wtaxpayers_2012)
tab _hh_ETNW_2012_bin_wide if $ifSRwide & $if_TVHS_cTNW & year==2012

capture drop _hh_ETNW_2012_bin_small
gen long  _hh_ETNW_2012_bin_small = 995 if inrange( hh_cTNW_2012-hh_N_wtaxpayers_2012*750000, -500000,-200000)
replace   _hh_ETNW_2012_bin_small = 998 if inrange( hh_cTNW_2012-hh_N_wtaxpayers_2012*750000, -200000,200000)
replace   _hh_ETNW_2012_bin_small = 1002 if inrange( hh_cTNW_2012-hh_N_wtaxpayers_2012*750000, 200000,999999999)

*-------------------------------------------------------------------------* 
* 			Event plot regressions
*-------------------------------------------------------------------------*

gen hh_logk3371 = log(1000*1.02^(year-2012) + hh_k3371)

gen byte hh_k3371_dummy = hh_k3371>0 if !missing(hh_k3371)

* Decompose giving participation into two terms, cumulaive _exit and cumulative _entry
gen byte _exit  = (hh_k3371==0)*(L.hh_k3371>0) if !missing(hh_k3371) & !missing(L.hh_k3371) & year>=2013
replace _exit = 0 if !missing(hh_k3371) & year==2012

gen byte _entry = (hh_k3371>0)*(L.hh_k3371==0) if !missing(hh_k3371) & !missing(L.hh_k3371) & year>=2013
replace _entry = 1 if !missing(hh_k3371) & year==2012 & hh_k3371>0

by LOPNR (year):  gen hh_k3371_cmlexit = sum(_exit) 	if year>=2012
by LOPNR (year):  gen hh_k3371_cmlentry = sum(_entry)	if year>=2012

* Sample split, ex-ante givers v. nongivers
gen hh_logk3371_pregiv 		= hh_logk3371 if hh_pre_giver == 1
gen hh_logk3371_prenongiv 	= hh_logk3371 if hh_pre_giver == 0

bysort hhid_2012 year: egen hh_rentinnt = total(rentinnt), mi
bysort LOPNR (year): egen double hh_rentinnt_2012 = total(hh_rentinnt*(year==2012)),mi
* subsample of liquid households (avg wtaxrate = 0.01, capitalize deposits at r=2.4%, avg intr rate in 2012 [ssb])
gen hh_logk3371_liquid = hh_logk3371 if 0.01*hh_MVH_secM_2012<0.10*hh_rentinnt_2012/0.024

* Different constants for log shifting
gen hh_logk3371_c100	= log(100 *1.02^(year-2012) 	+ hh_k3371 )
gen hh_logk3371_c500	= log(500* 1.02^(year-2012) 	+ hh_k3371 )
gen hh_logk3371_c1500	= log(1500*1.02^(year-2012) 	+ hh_k3371 )


gen log_g_pregiv = log_g if hh_pre_giver==1
gen byte g_pos_prenongiv = hh_gaver_tot_dummy if hh_pre_giver==0
gen byte g_pos = hh_gaver_tot_dummy
gen log_g_condpos = log_g if g_pos==1	

gen log_g_rel_c12 = log_g_rel if hh_gaver_rel_2012>0 & !missing(hh_gaver_rel_2012)
gen log_g_nonrel_c12 = log_g_nonrel if hh_gaver_nonrel_2012>0 & !missing(hh_gaver_nonrel_2012)
gen log_g_intl_c12 = log_g_intl if hh_gaver_intl_2012>0 & !missing(hh_gaver_intl_2012)
gen log_g_dom_c12 = log_g_dom if (hh_gaver_nonrel_2012-hh_gaver_intl_2012 > 0)

gen log_g_liquid = log_g if 0.011*hh_MVH_secM_2012*1000000<0.1*hh_rentinnt_2012/0.024  /* Max wtax effect is less than  2012 intr income*/

* sample split by giving frequency
gen log_g_regular 		= log_g if inlist(_N_posgiving_pre,6,7)
gen log_g_occasional 		= log_g if inlist(_N_posgiving_pre,2,3)
gen hh_logk3371_regular 	= hh_logk3371 if inlist(_N_posgiving_pre,6,7)
gen hh_logk3371_occasional	= hh_logk3371 if inlist(_N_posgiving_pre,1,2,3)


*=================================================================================*
* Summary Statistics
*=================================================================================*


*---------------------------------------*
* Main summary statistics table
*---------------------------------------*

* Focus on if_TVHS_TNW sample
preserve
keep if  ${if_TVHS_cTNW_2} & !missing(gaver_tot_w) & !missing(hh_wtax_b12) & !missing(hh_MVH_secM_2012) & !missing(age) & !missing(hh_loggrossinc) & hh_head_2012==1

gen hh_wtax_b12_condMVHS = hh_wtax_b12 if hh_MVH_secM_2012>0

eststo sumstats1a: estpost tabstat hh_gaver_tot_dummy  					, statistics(N mean 		)  columns(statistics)
eststo sumstats1b: estpost tabstat hh_gaver_tot_w 		if gaver_tot_w>0  	, statistics(N mean p25 p50 p75 )  columns(statistics)
eststo sumstats2a: estpost tabstat hh_wtax_dummy					, statistics(N mean 		)  columns(statistics)
eststo sumstats2b: estpost tabstat hh_wtax_b12			if hh_wtax_b12>0  	, statistics(N mean p25 p50 p75 )  columns(statistics)
eststo sumstats2bb:estpost tabstat hh_wtax_b12_condMVHS				 	, statistics(N mean p25 p50 p75 )  columns(statistics)
eststo sumstats2c: estpost tabstat hh_MWTR hh_AWTR		if hh_wtax_b12>0	, statistics(N mean 		)  columns(statistics)
eststo sumstats3:  estpost tabstat hh_MVH_priseccabM_2012  	if  year==2012		, statistics(N mean p25 p50 p75 )  columns(statistics)
eststo sumstats4:  estpost tabstat hh_MVH_secM_2012 		if hh_MVH_secM_2012>0 &  year==2012, statistics(N mean p25 p50 p75 )  columns(statistics)

eststo sumstats4b:  estpost tabstat hh_MNW 			if year==2012		, statistics(N mean p25 p50 p75 )  columns(statistics)
eststo sumstats5:  estpost tabstat age 				if year==2012		, statistics(N mean p25 p50 p75 )  columns(statistics)
eststo sumstats6:  estpost tabstat hh_grossinc 			if year==2012 		, statistics(N mean p25 p50 p75 )  columns(statistics)
eststo sumstats7:  estpost tabstat hh_N_wtaxpayers 		if year==2012		, statistics(N mean p25 p50 p75 )  columns(statistics)
restore

*[TABLE A.1]
foreach X in "1a" "1b" "2a" "2b" "2bb" "2c" "3" "4" "4b" "5" "6" "7" {
	local REPLACEORAPPEND = "append"
	if "`X'"=="1a" {
		local REPLACEORAPPEND = "replace"
	}
	
	local decimals = 2
	if inlist("`X'","1b","2b","2bb","4b","5","6") {
		local decimals = 0
	}
	if inlist("`X'","2c") {
		local decimals = 4
	}
	local cells = "count(fmt(%12.0fc)) Mean(fmt(%12.`decimals'fc)) p25(fmt(%12.`decimals'fc)) p50(fmt(%12.`decimals'fc)) p75(fmt(%12.`decimals'fc))"
	#delimit;	
	esttab sumstats`X',			
			cells("`cells'") 
			substitute(	hh_gaver_tot_dummy 	"1[Giving$ _{h,t} >0$]"
					hh_gaver_tot_w		"Giving$ _{h,t}$ if $ >0$"
					hh_wtax_dummy   	"$ 1[wtax_{h,t}>0]$"
					
					hh_wtax_b12_condMVHS    "$ wtax_{h,t}$ if MVHS$ _{h,2012}>0$"
					hh_wtax_b12		"$ wtax_{h,t}$ if $ >0$"
					hh_MVH_priseccabM_2012 	"MVH$ _{h,2012}$, MNOK"
					hh_MVH_secM_2012	"MVHS$ _{h,2012}$ if $ >0$, MNOK"
					age			"Age$ _{i,2012}$"
					hh_grossinc		"Gross income$ _{h,2012}$"
					hh_MWTR			"MWTR$ _{h,t}$"
					hh_AWTR			"AWTR$ _{h,t}$"
					hh_N_wtaxpayers		"Number of adults$ _h$"
					hh_MNW 			"Net wealth$ _{h,t}$"
					)	
			booktabs nolines gaps nonum nomtitle collabels(none) noobs frag `REPLACEORAPPEND'
			posthead()
			prefoot(); 
	#delimit cr	
}

tab _N_posgiving_pre if year ==2012 & ${if_TVHS_cTNW_2} & !missing(gaver_tot_w) & !missing(hh_wtax_b12) & !missing(hh_MVH_secM_2012) & !missing(age) & !missing(hh_loggrossinc) & hh_head_2012==1

*[FIGURE A.17]
#delimit;
hist _N_posgiving_pre 
	if year ==2012 & ${if_TVHS_cTNW_2} 
	& !missing(gaver_tot_w) & !missing(hh_wtax_b12) & !missing(hh_MVH_secM_2012) & !missing(age) & !missing(hh_loggrossinc) & hh_head_2012==1
	, freq discrete graphregion(color(white)) color(blue%80) lcolor(blue%60)
	  barwidth(0.75)
	  xtitle("Number of years with positive giving, 2006 to 2012")
	  ytitle("Number of households")
	  ylabel(0(100000)400000, format(%12.0g) glcolor(gs12%0) gmin gmax)
	  ymtick(0(50000)410000)
	  xlabel(0(1)6);
	  

#delimit cr

graph save ${figures}/wdon1_histogram_N_posgiving_pre.gph, replace
graph export ${figures}/wdon1_histogram_N_posgiving_pre.pdf,  as(pdf) replace

*-----------------------------------------*
* Summary statistics by whether MVHS>0
*-----------------------------------------*
* Focus on if_TVHS_TNW sample
preserve
keep if  ${if_TVHS_cTNW_2} & !missing(gaver_tot_w) & !missing(hh_wtax_b12) & !missing(hh_MVH_secM_2012) & !missing(age) & !missing(hh_loggrossinc) & year==2012 & hh_head_2012==1
gen hh_MNW_MNOK = hh_MNW/1000000
gen hh_grossinc_MNOK = hh_grossinc/1000000

foreach SUBSET of numlist 0 1 {
	if `SUBSET'==0 {
		global ifstatement = "hh_MVH_secM_2012 ==0"
	}
	if `SUBSET'==1 {
		global ifstatement = "hh_MVH_secM_2012 >0"
	}
	eststo sumstats1_`SUBSET': estpost tabstat hh_gaver_tot_dummy  		if $ifstatement		, statistics(N mean 		)  columns(statistics)
	eststo sumstats2_`SUBSET': estpost tabstat hh_gaver_tot_w 		if $ifstatement & gaver_tot_w>0  	, statistics(N mean p25 p50 p75 )  columns(statistics)
	eststo sumstats3_`SUBSET': estpost tabstat hh_wtax_dummy		if $ifstatement		, statistics(N mean 		)  columns(statistics)
	eststo sumstats4_`SUBSET': estpost tabstat hh_wtax_b12			if $ifstatement & hh_wtax_b12>0  	, statistics(N mean p25 p50 p75 )  columns(statistics)
	eststo sumstats5_`SUBSET': estpost tabstat hh_MVH_priseccabM_2012  	if $ifstatement 	, statistics(N mean p25 p50 p75 )  columns(statistics)
	eststo sumstats6_`SUBSET': estpost tabstat hh_MVH_secM_2012 		if $ifstatement		, statistics(N mean p25 p50 p75 )  columns(statistics)
	eststo sumstats7_`SUBSET': estpost tabstat hh_MNW_MNOK 			if $ifstatement		, statistics(N mean p25 p50 p75 )  columns(statistics)
	eststo sumstats8_`SUBSET': estpost tabstat age 				if $ifstatement		, statistics(N mean p25 p50 p75 )  columns(statistics)
	eststo sumstats9_`SUBSET': estpost tabstat hh_grossinc_MNOK 		if $ifstatement		, statistics(N mean p25 p50 p75 )  columns(statistics)
	eststo sumstats10_`SUBSET':estpost tabstat hh_N_wtaxpayers 		if $ifstatement		, statistics(N mean p25 p50 p75 )  columns(statistics)
}
restore

*[TABLE A.2]
foreach i of numlist 1(1)10 {
	#delimit;
	if `i' == 1  {; local decimals = 2; };
	if `i' == 2  {; local decimals = 0; };
	if `i' == 3  {; local decimals = 2; };
	if `i' == 4  {; local decimals = 0; };
	if `i' == 5  {; local decimals = 2; };
	if `i' == 6  {; local decimals = 2; };
	if `i' == 7  {; local decimals = 2; };
	if `i' == 8  {; local decimals = 0; };
	if `i' == 9  {; local decimals = 2; };
	if `i' == 10 {; local decimals = 2; };
	#delimit cr
	* count(fmt(%12.0fc)) dropped
	local cells = "Mean(fmt(%12.`decimals'fc)) p25(fmt(%12.`decimals'fc)) p50(fmt(%12.`decimals'fc)) p75(fmt(%12.`decimals'fc))"
	#delimit; 
	esttab  sumstats`i'_0 sumstats`i'_1,			
			cells("`cells'") 
			substitute(	hh_gaver_tot_dummy 	"1[Giving$ _{h,2012} >0$]"
					hh_gaver_tot_w		"Giving$ _{h,2012}$ if $ >0$"
					hh_wtax_dummy   	"$ 1[wtax_{h,2012}>0]$"
					
					hh_wtax_b12_condMVHS    "$ wtax_{h,t}$ if MVHS$ _{h,2012}>0$"
					hh_wtax_b12		"$ wtax_{h,t}$ if $ >0$"
					hh_MVH_priseccabM_2012 	"MVH$ _{h,2012}$, MNOK"
					hh_MVH_secM_2012	"MVHS$ _{h,2012}$, MNOK"
					age			"Age$ _{i,2012}$"
					hh_MNW_MNOK			"Net Wealth$ _{h,2012}$, MNOK"
					hh_grossinc_MNOK		"Gross income$ _{h,2012}$, MNOK)"
					hh_MWTR			"MWTR$ _{h,2012}$"
					hh_AWTR			"AWTR$ _{h,2012}$"
					hh_N_wtaxpayers		"Number of adults$ _h$"
					)	
			booktabs nolines nonum nomtitle collabels(none) noobs frag 
			posthead()
			prefoot(); 
	#delimit cr
}

gen loginterest_minus_logwtax = log(1+hh_rentinnt)-log(hh_wtax)
sum loginterest_minus_logwta if hh_gaver_tot>0 & hh_MVH_cabsecM_2012>0 & ${if_TVHS_cTNW_2} & !missing(gaver_tot_w) & !missing(hh_wtax_b12) & !missing(hh_MVH_secM_2012) & !missing(age) & !missing(hh_loggrossinc) & year>2012 & hh_head_2012==1,d
*-----------------------------------------------------------------------*
* Summary statistics by whether pre-period giving ... among MVHS>0
*-----------------------------------------------------------------------*
preserve
keep if  ${if_TVHS_cTNW_2} & !missing(gaver_tot_w) & !missing(hh_wtax_b12) ///
	& !missing(hh_MVH_secM_2012) & !missing(age) & !missing(hh_loggrossinc) ///
	& year==2012 & hh_head_2012==1 ///
	& hh_MVH_secM_2012>0
	
gen hh_MNW_MNOK = hh_MNW/1000000
gen hh_grossinc_MNOK = hh_grossinc/1000000


foreach SUBSET of numlist 0 1 {
	if `SUBSET'==0 {
		global ifstatement = "hh_pre_giver ==0"
	}
	if `SUBSET'==1 {
		global ifstatement = "hh_pre_giver ==1"
	}
	eststo sumstats2_`SUBSET': estpost tabstat hh_gaver_tot_w 		if $ifstatement 	, statistics(N mean p25 p50 p75 )  columns(statistics)
	eststo sumstats3_`SUBSET': estpost tabstat hh_wtax_dummy		if $ifstatement		, statistics(N mean 		)  columns(statistics)
	eststo sumstats4_`SUBSET': estpost tabstat hh_wtax_b12			if $ifstatement & hh_wtax_b12>0  	, statistics(N mean p25 p50 p75 )  columns(statistics)
	eststo sumstats5_`SUBSET': estpost tabstat hh_MVH_priseccabM_2012  	if $ifstatement 	, statistics(N mean p25 p50 p75 )  columns(statistics)
	eststo sumstats6_`SUBSET': estpost tabstat hh_MVH_secM_2012 		if $ifstatement		, statistics(N mean p25 p50 p75 )  columns(statistics)
	eststo sumstats7_`SUBSET': estpost tabstat hh_MNW_MNOK 			if $ifstatement		, statistics(N mean p25 p50 p75 )  columns(statistics)
	eststo sumstats8_`SUBSET': estpost tabstat age 				if $ifstatement		, statistics(N mean p25 p50 p75 )  columns(statistics)
	eststo sumstats9_`SUBSET': estpost tabstat hh_grossinc_MNOK 		if $ifstatement		, statistics(N mean p25 p50 p75 )  columns(statistics)
	eststo sumstats10_`SUBSET':estpost tabstat hh_N_wtaxpayers 		if $ifstatement		, statistics(N mean p25 p50 p75 )  columns(statistics)
}
restore

*[TABLE A.3]
foreach i of numlist 1(1)10 {
	#delimit;
	if `i' == 1  {; local decimals = 2; };
	if `i' == 2  {; local decimals = 0; };
	if `i' == 3  {; local decimals = 2; };
	if `i' == 4  {; local decimals = 0; };
	if `i' == 5  {; local decimals = 2; };
	if `i' == 6  {; local decimals = 2; };
	if `i' == 7  {; local decimals = 2; };
	if `i' == 8  {; local decimals = 0; };
	if `i' == 9  {; local decimals = 2; };
	if `i' == 10 {; local decimals = 2; };
	#delimit cr
	* count(fmt(%12.0fc)) dropped
	local cells = "Mean(fmt(%12.`decimals'fc)) p25(fmt(%12.`decimals'fc)) p50(fmt(%12.`decimals'fc)) p75(fmt(%12.`decimals'fc))"
	#delimit; 
	esttab  sumstats`i'_0 sumstats`i'_1,			
			cells("`cells'") 
			substitute(	hh_gaver_tot_dummy 	"1[Giving$ _{h,2012} >0$]"
					hh_gaver_tot_w		"Giving$ _{h,2012}$ "
					hh_wtax_dummy   	"$ 1[wtax_{h,2012}>0]$"
					
					hh_wtax_b12_condMVHS    "$ wtax_{h,t}$ if MVHS$ _{h,2012}>0$"
					hh_wtax_b12		"$ wtax_{h,t}$ if $ >0$"
					hh_MVH_priseccabM_2012 	"MVH$ _{h,2012}$, MNOK"
					hh_MVH_secM_2012	"MVHS$ _{h,2012}$, MNOK"
					age			"Age$ _{i,2012}$"
					hh_MNW_MNOK			"Net Wealth$ _{h,2012}$, MNOK"
					hh_grossinc_MNOK		"Gross income$ _{h,2012}$, MNOK"
					hh_MWTR			"MWTR$ _{h,2012}$"
					hh_AWTR			"AWTR$ _{h,2012}$"
					hh_N_wtaxpayers		"Number of adults$ _h$"
					)	
			booktabs nolines nonum nomtitle collabels(none) noobs frag 
			posthead()
			prefoot(); 
	#delimit cr
}

sum hh_gaver_tot_dummy if year>2012 & hh_pre_giver==0 & ${if_TVHS_cTNW_2} & !missing(gaver_tot_w) & !missing(hh_wtax_b12) & !missing(hh_MVH_secM_2012) & !missing(age) & !missing(hh_loggrossinc)

*---------------------------------------*
* Misc. statistics for use in text
*---------------------------------------*
* share of secondary home owners paying wtax as of 2015 (mid sample)
sum hh_wtax_dummy if hh_head==1 & year==2015
sum hh_wtax_dummy if hh_head==1 & year==2015 & !missing(hh_TVH_pri) & hh_TVH_pri>0
sum hh_wtax_dummy if hh_head==1 & year==2015 & !missing(hh_TVH_sec) & hh_TVH_sec>0
sum hh_wtax_dummy if hh_head==1 & year==2018 & !missing(hh_TVH_sec) & hh_TVH_sec>0

* marginal share of secondary housing subj to wealth tax
gen _margshare = min(max(hh_TNW-hh_N_wtaxpayers*wtax_threshold,0),hh_TVH_sec)/hh_TVH_sec
replace _margshare = 0 if !missing(hh_TNW) & missing(_margshare)
sum _margshare if hh_head==1 & year==2018 & !missing(hh_TVH_sec) & hh_TVH_sec>0 & hh_wtax_dummy==1


*=============================================================================*
* Bunching analysis & plots
*=============================================================================*
capture drop gaver_intl
gen gaver_intl = gaver_aid + gaver_climate + gaver_humrig

foreach givingvar in  "gaver_tot"  "gaver_nonrel" "gaver_intl"  "gaver_rel" {
	global givingvar = "`givingvar'"
	capture drop EG
	if "$givingvar" == "gaver_rel" { 
		/* exclusively religious giving*/
		gen double EG = gaver_tot-don_cap if gaver_tot==gaver_rel
		global appendfilename = "_rel"
	}
	if "$givingvar" == "gaver_nonrel" { 
		/* exclusively non-religious giving*/
		gen double EG = gaver_tot-don_cap if gaver_rel==0
		global appendfilename = "_nonrel"
	}
	if "$givingvar" == "gaver_tot" {
		gen double EG = gaver_tot-don_cap
		global appendfilename = ""
	}
	if "$givingvar" == "gaver_intl" {
		gen double EG = gaver_tot-don_cap if gaver_intl==gaver_tot
		global appendfilename = "_intl"
	}
	if "$givingvar" == "gaver_dom" {
		gen double EG = gaver_tot-don_cap if gaver_dom==gaver_tot
		global appendfilename = "_dom"
	}

	*-----------------------------------------------------------------*
	* Basic variables for bunch_count
	*-----------------------------------------------------------------*

	capture drop *G*bin*
	gen long _EG_bin = floor(EG/100)*100
	bysort _EG_bin: gen _EG_bin_N = _N
	by _EG_bin: gen _EG_bin_n = _n

	*-----------------------------------------------------------------*
	* Some stuff to account for round-number bunching
	*-----------------------------------------------------------------*

	gen long _G_bin = floor(gaver_tot/100)*100 if !missing(EG)

	sort _G_bin
	foreach year of numlist 2012(1)2018 {
		by _G_bin: egen long _G_bin_`year'_N = total(year==`year')
	}

	bysort 	year _G_bin: gen long _G_bin_year_N = _N
	by 	year _G_bin: gen byte _G_bin_year_unique = _n==1

	* Assign a placebo # of observations
	gen long _G_bin_year_placeboN = .
	replace _G_bin_year_placeboN = (				_G_bin_2014_N + _G_bin_2015_N + _G_bin_2016_N + _G_bin_2017_N + _G_bin_2018_N)/5 if year == 2012
	replace _G_bin_year_placeboN = (				_G_bin_2014_N + _G_bin_2015_N + _G_bin_2016_N + _G_bin_2017_N + _G_bin_2018_N)/5 if year == 2013
	replace _G_bin_year_placeboN = (_G_bin_2012_N + _G_bin_2013_N + 		_G_bin_2015_N + _G_bin_2016_N + _G_bin_2017_N + _G_bin_2018_N)/6 if year == 2014
	replace _G_bin_year_placeboN = (_G_bin_2012_N + _G_bin_2013_N + _G_bin_2014_N + 		_G_bin_2016_N + _G_bin_2017_N + _G_bin_2018_N)/6 if year == 2015
	replace _G_bin_year_placeboN = (_G_bin_2012_N + _G_bin_2013_N + _G_bin_2014_N +	_G_bin_2015_N +			_G_bin_2017_N + _G_bin_2018_N)/6 if year == 2016
	replace _G_bin_year_placeboN = (_G_bin_2012_N + _G_bin_2013_N + _G_bin_2014_N +	_G_bin_2015_N + _G_bin_2016_N + 		_G_bin_2018_N)/6 if year == 2017
	replace _G_bin_year_placeboN = (_G_bin_2012_N + _G_bin_2013_N + _G_bin_2014_N +	_G_bin_2015_N + _G_bin_2016_N + _G_bin_2017_N 		     )/6 if year == 2018

	* For each bin, for each year; Also want to count # obs in neighboring bins
	bysort  _G_bin_year_unique year (_G_bin)  : gen _G_binneighbors_year_N = ceil(0.5*(_G_bin_year_N[_n-1] + _G_bin_year_N[_n+1]))

	* Now, or non-round-number bins, just use bin count
	gen byte _G_binisaroundnumb_year = mod(_G_bin,1000)==0
	gen int _G_bin_N_year_adj = _G_bin_year_N if _G_binisaroundnumb_year==0

	* But for round-number bins, adjust
	replace _G_bin_N_year_adj = _G_bin_year_N -_G_bin_year_placeboN + _G_binneighbors_year_N  if _G_binisaroundnumb_year==1
	replace _G_bin_N_year_adj = max(_G_binneighbors_year_N,_G_bin_N_year_adj)  if _G_binisaroundnumb_year==1 

	* Now count adjusted number of obs per bin across years
	bysort _EG_bin: egen _EG_bin_N_adj = total(_G_bin_N_year_adj * _G_bin_year_unique),mi 

	* Needs to stay sorted by _EG_bin for bunch_count

	*-------------------------------*

	*-------------------------------*
	* Needed to calculate relative excess mass / elasticity
	sum don_cap if inrange(_EG_bin,-4000,3900) & !missing(_EG_bin_N) 
	global mean_don_cap = `r(mean)'
	sum don_rate if inrange(_EG_bin,-4000,3900) & !missing(_EG_bin_N) 
	global mean_don_rate = `r(mean)'
	
	* Needed for text box placement (need y-axis max)
	sum _EG_bin_N_adj if inrange(_EG_bin,-4000,3900) & !missing(_EG_bin_N) 
	global Nmax = `r(max)'
	
	* Needed for text box placement (need y-axis max)
	sum _EG_bin_N_adj if _EG_bin_n==1 & inrange(_EG_bin,-4000,3900)
	global Nmax = `r(max)'		
	global Nmax_divby_2 = ${Nmax}/2

	*===> Adjusted
		preserve
		keep if _EG_bin_n==1
		capture bunch_count _EG_bin _EG_bin_N_adj if _EG_bin_n==1, ig_low(-40) ig_high(39) low_bunch(-1) high_bunch(1) binwidth(100) plot(1) nboot(200)
		restore
		graph save ${figures}/wdon1_bunch_count_pooled_adj${appendfilename}_RAW.gph, replace
		graph use ${figures}/wdon1_bunch_count_pooled_adj${appendfilename}_RAW.gph
		global b_est =subinstr(subinstr("abc"+string(round(b,0.01)),"abc.","abc0.",.),"abc","",.)
		global B_est = string(round(b*100/${mean_don_cap},0.01)) /* 100 is binwidth */
		global b_se = subinstr(subinstr("abc"+string(round(b_se,0.01)),"abc.","abc0.",.),"abc","",.)
		global Kbar = string(round(${mean_don_cap}))
		global tau_g =subinstr(subinstr("abc"+string(round(${mean_don_rate},0.0001)),"abc.","abc0.",.),"abc","",.)
		global e_est_ = (b*100/${mean_don_cap})/(-log(1-${mean_don_rate}))
		global e_est  ="-" + substr(subinstr(subinstr("abc"+string(round(${e_est_},0.01)),"abc.","abc0.",.),"abc","",.) + "0",1,4) /* keeps leading and lagging zeros..*/
		global e_se_ = (b_se*100/${mean_don_cap})/(-log(1-${mean_don_rate}))
		global e_se  =substr(subinstr(subinstr("abc"+string(round(${e_se_},0.01)),"abc.","abc0.",.),"abc","",.) + "0",1,4)
		
		* Fix 
		gr_edit .style.editstyle margin(medlarge) editcopy
		gr_edit .yaxis1.style.editstyle majorstyle( gridstyle( linestyle(color(gs14%0)) draw_min(yes) draw_max(yes))) editcopy
		gr_edit .plotregion1.plot1.style.editstyle marker(size(0.5)  symbol(circle) linestyle(color(blue)) fillcolor(blue)	) editcopy
		gr_edit .plotregion1.plot1.style.editstyle line(color(blue))	editcopy
		gr_edit .plotregion1.plot2.style.editstyle line(color(green) width(medthick)) editcopy
		gr_edit .plotregion1._xylines[1].style.editstyle linestyle(color(gs4%15) width(vthick)) editcopy
		gr_edit .title.text = {"{stSerif: Panel B: Bunching at Deduction Cap}"} 
		gr_edit .subtitle.text = {    "{stSerif:Estimated e=${e_est} (se=${e_se})}"  }
		*gr_edit .subtitle.text.Arrpush  "{stSerif:b=${b_est} (${b_se}), mean K{sub:t}=${Kbar}, {&tau}{sub:g,t}=${tau_g}}"  
		gr_edit .title.style.editstyle size(medlarge) color(black) margin(vsmall) editcopy
		gr_edit .subtitle.style.editstyle size(medium) color(blue) margin(small) editcopy
		
		gr_edit ..plotregion1.AddTextBox added_text editor  ${Nmax_divby_2} 20
		gr_edit ..plotregion1.added_text[1].style.editstyle linegap(0.5) color(gs5) editcopy
		gr_edit ..plotregion1.added_text[1].text = {}
		gr_edit ..plotregion1.added_text[1].text.Arrpush `"b=${b_est} (se=${b_se})"'
		gr_edit ..plotregion1.added_text[1].text.Arrpush `"mean K{sub:t}=${Kbar}"'
		gr_edit ..plotregion1.added_text[1].text.Arrpush `"mean {&tau}{sup:g}{sub:t}=${tau_g}"'


			foreach b of numlist 9(-1)1 {
				local oldnum = round(-40 +(`b'-1)*10)
				local newnum = round(`oldnum'*100)
				disp "`b' , `oldnum' , `newnum'"
				gr_edit .xaxis1.major.num_rule_ticks = 0
				gr_edit .xaxis1.edit_tick `b' `oldnum' "`newnum'", tickset(major)
			}

		gr_edit .yaxis1.title.text = {"Adjusted Frequency"}

		gr_edit .xaxis1.title.text = {"Charitable Giving in Excess of Cap (NOK)"}

		gr_edit .style.editstyle boxstyle(linestyle(color(white))) editcopy
		
		graph save ${figures}/wdon1_bunch_count_pooled_adj${appendfilename}.gph, replace
}



* Export individual plots as pdf without titles
foreach X in "" "_rel" "_nonrel" "_intl"   {
	global X = "`X'"
	graph use ${figures}/wdon1_bunch_count_pooled_adj${X}.gph
	gr_edit .title.text = {} 
	gr_edit .plotregion1.plot1.style.editstyle line(width(medthick)) editcopy
	gr_edit .plotregion1.plot2.style.editstyle line(width(thick)) editcopy
	graph export ${figures}/wdon1_bunch_count_pooled_adj${X}.pdf, as(pdf) replace
}


*[FIGURE 5]

gr combine ${figures}/wdon1_bunch_count_pooled_adj_rel.gph ${figures}/wdon1_bunch_count_pooled_adj_nonrel.gph, ///
	rows(1) graphregion(color(white) margin(tiny tiny tiny tiny)) xsize(2) ysize(1)   imargin(1 4 1 1)
gr_edit .plotregion1.graph1.yaxis1.reset_rule 0 50000 10000 , tickset(major) ruletype(range)	
gr_edit .plotregion1.graph2.yaxis1.reset_rule 0 5000 1000 , tickset(major) ruletype(range)
gr_edit .plotregion1.graph2.style.editstyle boxstyle(linestyle(color(white))) editcopy
gr_edit .plotregion1.graph2.title.text = {"{stSerif: Panel B: Nonreligious Giving}"} 
gr_edit .plotregion1.graph1.title.text = {"{stSerif: Panel A: Religious Giving}"} 
	gr_edit .plotregion1.graph1.plotregion1.plot1.style.editstyle line(width(medthick)) editcopy
	gr_edit .plotregion1.graph1.plotregion1.plot2.style.editstyle line(width(thick)) editcopy
	gr_edit .plotregion1.graph2.plotregion1.plot1.style.editstyle line(width(medthick)) editcopy
	gr_edit .plotregion1.graph2.plotregion1.plot2.style.editstyle line(width(thick)) editcopy
gr_edit plotregion1.graph1.title.style.editstyle size(large) color(black) margin(vsmall) editcopy
gr_edit plotregion1.graph2.title.style.editstyle size(large) color(black) margin(vsmall) editcopy
* may want to manually shift position of textbox in graph2 here
graph export ${figures}/wdon1_bunch_count_pooled_adj_comb_rel_nonrel.pdf, as(pdf) replace

*-------------------------------------------------*
* [FIGURE 4] Bunching at the tax-deduction cap for charitable giving
*-------------------------------------------------*
#delimit;
hist gaver_tot if inrange(gaver_tot,8000,55000), 
	freq width(200) graphregion(color(white)) color(blue) 
	xtitle("Charitable Giving (NOK)")
	xlabel(10000(10000)55000) xmtick(5000(5000)55000)
	title("{stSerif: Panel A: Pooled Histogram}")
	ylabel(0(10000)40000, format(%12.0g) glcolor(gs12%0) gmin gmax) 
	ymtick(5000(5000)50000)
	xline(12000, lcolor(gray%20) lwidth(vvthick))
	xline(16800, lcolor(gray%20) lwidth(vvthick))
	xline(20000, lcolor(gray%20) lwidth(vvthick))
	xline(25000, lcolor(gray%20) lwidth(vvthick))
	xline(30000, lcolor(gray%20) lwidth(vvthick))
	xline(40000, lcolor(gray%20) lwidth(vvthick))
	xline(50000, lcolor(gray%20) lwidth(vvthick));
#delimit cr
gr_edit .style.editstyle margin(medlarge) editcopy
gr_edit .subtitle.text = {" "}
gr_edit .title.style.editstyle size(medlarge) color(black) margin(vsmall) editcopy	
	
graph save ${figures}/wdon1_histogram_pooled.gph, replace
graph export ${figures}/wdon1_histogram_pooled.pdf,  as(pdf) replace



**** Combined w/ adjusted
gr combine ${figures}/wdon1_histogram_pooled.gph ${figures}/wdon1_bunch_count_pooled_adj.gph, ///
	rows(1) graphregion(color(white) margin(tiny tiny tiny tiny)) xsize(2) ysize(1) ycommon  imargin(1 1 1 1)
gr_edit .plotregion1.graph2.style.editstyle boxstyle(linestyle(color(white))) editcopy
gr_edit .plotregion1.graph2.plotregion1.plot1.style.editstyle line(width(medthick)) editcopy
gr_edit .plotregion1.graph2.plotregion1.plot2.style.editstyle line(width(thick)) editcopy
gr_edit plotregion1.graph1.title.style.editstyle size(vlarge) color(black) margin(vsmall) editcopy
gr_edit plotregion1.graph2.title.style.editstyle size(vlarge) color(black) margin(vsmall) editcopy


graph export ${figures}/wdon1_overview_donation_bunching.pdf, as(pdf) replace



*-------------------------------------------------*
* Elasticity Heterogeneity Regs (bunching)
*-------------------------------------------------*

gen byte Ind = inrange(gaver_tot-don_cap,-100,99) 
gen int Ind100 = Ind*100 /* To get coefs in % */

gen _g1 = gaver_tot
gen _g2 = gaver_tot^2
gen _g3 = gaver_tot^3


gen _hh_wtax_b12_MNOK = hh_wtax_b12/1000000


foreach Yvar of varlist Ind  {
	global Yvar = "`Yvar'"
	eststo reg${Yvar}_1: reghdfe ${Yvar}100 alder 	  	 		  			if $ifSRwide & $if_RKD_TNW , absorb(hh_N_wtaxpayers year##c.(_g1 _g2 _g3) )  cluster(hhid)
	eststo reg${Yvar}_2: reghdfe ${Yvar}100 	hh_loggrossinc   				if $ifSRwide & $if_RKD_TNW , absorb(hh_N_wtaxpayers year##c.(_g1 _g2 _g3) )  cluster(hhid)
	eststo reg${Yvar}_3: reghdfe ${Yvar}100 			hh_logMNW			if $ifSRwide & $if_RKD_TNW , absorb(hh_N_wtaxpayers year##c.(_g1 _g2 _g3) )  cluster(hhid)
	eststo reg${Yvar}_4: reghdfe ${Yvar}100 alder 	hh_loggrossinc  hh_logMNW  			if $ifSRwide & $if_RKD_TNW , absorb(hh_N_wtaxpayers year##c.(_g1 _g2 _g3) )  cluster(hhid)
}


*[TABLE 2]
#delimit;
esttab  regInd_1 regInd_2 regInd_3 regInd_4 
	, 
	b(4) se(4)
	varwidth(30) stats(FS_F r2 N, fmt(2 4 0 )) 
	 
	star(* 0.05 ** 0.01 *** 0.001)  nocons
	substitute ( 
			alder "Age"
			hh_loggrossinc "log(Household Gross Income)"
			hh_wtax_dummy "1[wealth tax$ >$0]"
			_hh_wtax_b12_MNOK "wealth tax (MNOK)"
			hh_logMNW "log(Market-value Net Wealth)"
			r2 "R $^2$"
	)
	booktabs nolines nonum nomtitle nogaps frag replace;
#delimit cr

* also want unconditional prob
sum Ind if $ifSRwide & $if_RKD_TNW
sum Ind if $ifSRwide & $if_RKD_TNW & G_share_rel>0.999
sum Ind if $ifSRwide & $if_RKD_TNW & G_share_rel<0.001



*========================================================================*
* Behavior Around Wealth Tax Threshold (RKD)
*========================================================================*

capture drop _hh_ETNW_bin* _dontomit
gen long _hh_ETNW_bin = floor(hh_cETNW/50000)*50000/1000
gen long _hh_ETNW_bin_sh = 2000 + _hh_ETNW_bin if inrange(_hh_ETNW_bin,-1000,1000-1)
gen _dontomit = _hh_ETNW_bin!=-250
tab _hh_ETNW_bin_sh _dontomit

capture drop _hh_ETNW_pos
gen _hh_ETNW_pos = hh_ETNW>0 if !missing(hh_ETNW)


global s=1
global Yvar = "hh_gaver_tot_w"
global YvarRegName = subinstr("$Yvar",".","_",.)

reghdfe ${Yvar} _hh_ETNW_bin_sh#c._dontomit alder hh_loggrossinc  if $ifSRwide & $if_RKD_TNW & hh_head==1, ///
	absorb(hh_N_wtaxpayers year )  cluster(hhid)
	
		mat regtable =r(table)'
		mat r${s}_${YvarRegName} = regtable[1..rownumb(regtable,"2950._hh_ETNW_bin_sh#c._dontomit"),1..6]
		
		
		capture drop r${s}_${YvarRegName}*
		svmat2  r${s}_${YvarRegName}, rnames(r${s}_${YvarRegName}_rn) 
		replace r${s}_${YvarRegName}_rn = substr(r${s}_${YvarRegName}_rn,1,4)
		replace r${s}_${YvarRegName}_rn = subinstr(r${s}_${YvarRegName}_rn,".","",.)
		replace r${s}_${YvarRegName}_rn = subinstr(r${s}_${YvarRegName}_rn,"_","",.)
		replace r${s}_${YvarRegName}_rn = subinstr(r${s}_${YvarRegName}_rn,"b","",.)
		replace r${s}_${YvarRegName}_rn = subinstr(r${s}_${YvarRegName}_rn,"o","",.)
		gen long r${s}_${YvarRegName}_xb = (real(r${s}_${YvarRegName}_rn)-2000)
		drop r${s}_${YvarRegName}2 r${s}_${YvarRegName}3  r${s}_${YvarRegName}4 
		

reghdfe ${Yvar} c.hh_ETNW c._hh_ETNW_pos c.hh_ETNW#c._hh_ETNW_pos  alder hh_loggrossinc 	///
	if !missing(_hh_ETNW_bin_sh) & $ifSRwide & $if_RKD_TNW & hh_head==1, ///
	absorb(hh_N_wtaxpayers year )  cluster(hhid)
	
	mat s${s}_${YvarRegName}= r(table)'
	
	* y = a + b*x, when x <0
	* Use this to fhistind intercept for graph: a = E[scatter points] - b*E[x]
	global b = s${s}_${YvarRegName}[1,1]
	global c = s${s}_${YvarRegName}[3,1]
	global d = s${s}_${YvarRegName}[2,1]
	#delimit;
	sum r${s}_${YvarRegName}1   if r${s}_${YvarRegName}_xb<0; local mean_y_left = `r(mean)';
	sum r${s}_${YvarRegName}_xb if r${s}_${YvarRegName}_xb<0; local mean_x_left = `r(mean)';
	#delimit cr
	
	global a = `mean_y_left' - ${b}*`mean_x_left'
	
	
	global dslope_y 	= string(round(${c}*1000*1000,0.1))
	global dslope_y_se 	= string(round(s${s}_${YvarRegName}[3,2]*1000*1000,0.1))
	global disc_y 		= string(round(${d},0.1))
	global disc_y_se 	= string(round(s${s}_${YvarRegName}[2,2],0.1))

* Figure with fitted lines
#delimit;
twoway  (scatter r${s}_${YvarRegName}1 r${s}_${YvarRegName}_xb, msize(small) mcolor(green) yaxis(1))
	(function y = ${a} + ${b}*x*1000       	 	, range(-1000 0   ) lcolor(green) lpattern(dash) yaxis(1))
	(function y = ${a} + ${b}*x*1000	 	 , range(  0 1000 ) lcolor(gs12) lpattern(dash) yaxis(1))
	(function y = ${a} + (${b}+${c})*x*1000+ ${d}	, range(  0 1000 ) lcolor(green) lpattern(dash) yaxis(1))
	, graphregion(color(white)) 
	ylabel(,glcolor(gs14) glwidth(vthin) gmin gmax)
	ylabel(-400(200)400)
	ytitle("Charitable Giving (NOK)", axis(1) color(black))
	legend(off)
	subtitle("Discontinuity = ${disc_y} (${disc_y_se}), dSlope=${dslope_y} (${dslope_y_se}) in MNOK")
	xtitle("Taxable Net Wealth - Threshold (NOK 1,000)")
	xline(0, lcolor(gs8) lpattern(dashed));
#delimit cr

graph save ${figures}/wdon1_RKD_${YvarRegName}.gph, replace
graph export ${figures}/wdon1_RKD_${YvarRegName}.pdf, replace as(pdf)

* Figure WITHOUT fitted lines
local endash = ustrunescape("\u2013")
#delimit;
twoway  (scatter r${s}_${YvarRegName}1 r${s}_${YvarRegName}_xb, msize(small) mcolor(blue) msymbol(Oh) yaxis(1))
	
	, graphregion(color(white)) 
	ylabel(,glcolor(gs14) glwidth(vthin) gmin gmax)
	ylabel(-400(200)400)
	ytitle("Charitable Giving (NOK)", axis(1) color(black))
	legend(off)
	xtitle("Taxable Net Wealth `endash' Threshold (NOK 1,000)")
	xline(0, lcolor(gs8) lpattern(dashed))
	xsize(1.75) ysize(1);
#delimit cr

graph save ${figures}/wdon1_RKD_${YvarRegName}_nofitted.gph, replace
graph export ${figures}/wdon1_RKD_${YvarRegName}_nofitted.pdf, replace as(pdf)

*===================================================================*
* Reduced-form regressions, event plots
*===================================================================*


foreach Yvar of varlist hh_logk3371_occasional   hh_logk3371_regular  hh_gaver_tot  hh_logk3371_c100 hh_logk3371_c500 hh_logk3371_c1500 hh_logk3371  hh_logk3371_liquid  hh_logk3371_pregiv hh_logk3371_prenongiv  hh_k3371_dummy hh_k3371_cmlexit hh_k3371_cmlentry    hh_AWTR hh_MWTR hh_logk3371 hh_wtax_b12 hh_k3371_w hh_gaver_tot_w {
	global Yvar = "`Yvar'"
	
	global year_start = 2010
	global year_end = 2018
	global base_control = "c.hh_MVH_priseccabM_2012"
	global inc_control = "hh_loggrossinc_2012"
	if inlist("${Yvar}","hh_k3371_w","hh_logk3371","hh_k3371_dummy","hh_k3371_entry","hh_k3371_exit","hh_k3371_cmlentry","hh_k3371_cmlexit","hh_logk3371_pregiv","hh_logk3371_prenongiv") {
		global year_start=2006
	}
	if inlist("${Yvar}","hh_logk3371_regular","hh_logk3371_occasional") {
		global year_start=2012
	}
	if inlist("${Yvar}","hh_logk3371_c100", "hh_logk3371_c500", "hh_logk3371_c2000") {
		global year_start=2006
	}
	if inlist("${Yvar}","hh_logk3371_liquid") {
		global year_start=2006
	}

	if inlist("${Yvar}","hh_wtax_b12","hh_MWTR","hh_AWTR") {
		global year_start=2010
	}
	
	disp "$year_start"

	#delimit;
	reghdfe	$Yvar
			year#c.hh_MVH_secM_2012#c.not2012
			c.hh_MVH_secM_2012 
			
		if $if_TVHS_cTNW_2 & hh_head_2012==1
		& inrange(year,${year_start},${year_end})
		, cluster(hhid_2012)
		absorb(year##c.(	$base_control	
					hh_TNW_2012	hh_TNW_2012_sq 	hh_TNW_2012_cu
					age_2012 age2_2012 age3_2012 famsize_2012 famsize2_2012
					${inc_control}
					hh_N_wtaxpayers_2012
					hh_MVH_cabsec_dummy_2012
				)
			)
		;
	#delimit cr
	estimates save ${temp_estimates}/event_${Yvar}, replace
	mat T${Yvar} = r(table)'
}



mat Shh_wtax_b12 = Thh_wtax_b12[1..9,1..6]
mat Shh_AWTR = Thh_AWTR[1..9,1..6]
mat Shh_k3371_w = Thh_k3371_w[1..13,1..6]
mat Shh_logk3371 = Thh_logk3371[1..13,1..6]
mat Shh_gaver_tot_w = Thh_gaver_tot_w[1..7,1..6]
mat Shh_gaver_tot = Thh_gaver_tot[1..7,1..6]

mat Shh_logk3371_regular = Thh_logk3371_regular[1..7,1..6]
mat Shh_logk3371_occasional = Thh_logk3371_occasional[1..7,1..6]

mat Shh_k3371_cmlentry= Thh_k3371_cmlentry[1..7,1..6]
mat Shh_k3371_cmlexit= Thh_k3371_cmlexit[1..7,1..6]
mat Shh_k3371_dummy= Thh_k3371_dummy[1..13,1..6]

mat Shh_logk3371_pregiv= Thh_logk3371_pregiv[1..13,1..6]
mat Shh_logk3371_prenongiv= Thh_logk3371_prenongiv[1..7,1..6]
mat Shh_logk3371_liquid= Thh_logk3371_liquid[1..7,1..6]

foreach Yvar of varlist  hh_logk3371_occasional  hh_logk3371_regular hh_gaver_tot_w  hh_gaver_tot hh_k3371_dummy hh_MWTR hh_logk3371_pregiv hh_logk3371_prenongiv  hh_k3371_cmlexit hh_k3371_cmlentry   hh_logk3371 hh_AWTR  hh_wtax_b12 hh_k3371_w  {
	
	
	capture drop S`Yvar'*
	svmat2  S`Yvar', rnames(S`Yvar'_rn) 
		
		replace S`Yvar'_rn = substr(S`Yvar'_rn,1,4)
		replace S`Yvar'_rn = subinstr(S`Yvar'_rn,".","",.)
		replace S`Yvar'_rn = subinstr(S`Yvar'_rn,"_","",.)
		replace S`Yvar'_rn = subinstr(S`Yvar'_rn,"b","",.)
		replace S`Yvar'_rn = subinstr(S`Yvar'_rn,"o","",.)
		gen long S`Yvar'_xb = (real(S`Yvar'_rn))
		drop S`Yvar'2 S`Yvar'3  S`Yvar'4 
}

capture drop x x_shifted
gen x = Shh_k3371_w_xb
gen x_shifted = Shh_k3371_w_xb + 0.15

gen Shh_k3371_w_xb_shifted = Shh_k3371_w_xb - 0.15
gen Shh_logk3371_xb_shifted = Shh_logk3371_xb - 0.15

*[FIGURE A.7] First-stage and Reduced-form Plot 
#delimit;
twoway  
	(connected Shh_k3371_w1 		Shh_k3371_w_xb_shifted if Shh_k3371_w_xb_shifted>=2006.5,		color(blue)    yaxis(1) msymbol(Sh))
	(rcap Shh_k3371_w6 Shh_k3371_w5 	Shh_k3371_w_xb_shifted if Shh_k3371_w_xb_shifted>=2006.5,		color(blue%50) yaxis(1))
	(connected Shh_wtax_b121 		Shh_wtax_b12_xb,		color(orange) 	yaxis(2) msymbol(Sh) lpattern(dash) )
	(rcap Shh_wtax_b126 Shh_wtax_b125 	Shh_wtax_b12_xb,		color(orange%50) yaxis(2) lpattern(dash) )
	
	, graphregion(color(white))
	ylabel(-3000(1000)3000, axis(2) nogrid)
	ylabel(-60(20)60, axis(1) gmin gmax glcolor(gs14))
	
	xlabel(200(2)2018)
	xline(2012, lcolor(gs4%15) lwidth(vvvthick))
	ytitle("Wealth Tax Amount (NOK, 2012=0)", axis(2) color(orange))
	ytitle("Amount of Charitable Giving (NOK, 2012=0)", axis(1) color(blue))
	legend(order(1 "Reduced-form Effect on Giving" 3 "First-stage Effect on Wealth Taxes" ) ring(0) pos(11) region(lcolor(gs12)) margin(r=8) bmargin(1 0.5 0.5 4) cols(1))
	xtitle("Year")
	xsize(1.75) ysize(1);
#delimit cr
graph save ${figures}/wdon1_TVHS_IV_eventplot_RF_FS.gph, replace
graph export ${figures}/wdon1_TVHS_IV_eventplot_RF_FS.pdf, as(pdf) replace

*[FIGURE 2] Reduced-form Plot log(Giving)
replace Shh_logk33715 = Shh_logk33711 if Shh_logk3371_xb==2012
replace Shh_logk33716 = Shh_logk33711 if Shh_logk3371_xb==2012
#delimit;
twoway  
	(connected Shh_logk33711 		Shh_logk3371_xb if Shh_logk3371_xb_shifted>=2005.5,		lcolor(blue)    yaxis(1) msymbol(s) mcolor(blue%75) lwidth(medthick))
	(line Shh_logk33716 	Shh_logk3371_xb if Shh_logk3371_xb_shifted>=2005.5,		color(blue%50) yaxis(1) lpattern(longdash))
	(line Shh_logk33715 	Shh_logk3371_xb if Shh_logk3371_xb_shifted>=2005.5,		color(blue%50) yaxis(1) lpattern(longdash))
	ylabel(-0.01(0.005)0.005, axis(1) gmin gmax glcolor(gs14))
	xlabel(2007(2)2018)
	xmlabel(2007(1)2018, nolabel)
	xline(2012, lcolor(gs4%15) lwidth(vvvthick))
	ytitle("log of Charitable Giving (2012=0)", axis(1) color(blue))
	legend(off)
	xtitle("Year")
	xsize(1.75) ysize(1);
#delimit cr
graph save ${figures}/wdon1_TVHS_IV_eventplot_RF_log_g.gph, replace
graph export ${figures}/wdon1_TVHS_IV_eventplot_RF_log_g.pdf, as(pdf) replace


* [FIGURE A.8] Compare k3371 with totals from giving data
gen Shh_gaver_tot_w_xb_shifted = Shh_gaver_tot_w_xb - 0.15

#delimit;
twoway  
	(connected Shh_k3371_w1 		Shh_k3371_w_xb_shifted if Shh_k3371_w_xb_shifted>=2007,		color(blue)    yaxis(1) msymbol(Sh))
	(rcap Shh_k3371_w6 Shh_k3371_w5 	Shh_k3371_w_xb_shifted if Shh_k3371_w_xb_shifted>=2007,		color(blue%50) yaxis(1))
	(scatter Shh_gaver_tot_w1 		Shh_gaver_tot_w_xb_shifted ,	color(red)    yaxis(1) msymbol(Oh) lpattern(shortdash))
	(rcap Shh_gaver_tot_w6 Shh_gaver_tot_w5 Shh_gaver_tot_w_xb_shifted ,	color(red%50) yaxis(1))
	
	, graphregion(color(white))
	/*ymlabel(-0.01(0.0025)0.01, grid glcolor(gs14) glwidth(vthin) gmin gmax nolabel)*/
	ylabel(-60(20)60, axis(1) gmin gmax glcolor(gs14))
	
	xlabel(2007(2)2018)
	xline(2012, lcolor(gs4%15) lwidth(vvvthick))
	ytitle("Amount of Charitable Giving (NOK, 2012=0)", axis(1) color(blue))
	legend(order(1 "Giving (long panel)" 3  "Giving (short panel)"  ) ring(0) pos(11) region(lcolor(gs12)) margin(r=8) bmargin(1 0.5 0.5 4) cols(1))
	xtitle("Year")
	xsize(1.75) ysize(1);
#delimit cr
graph save ${figures}/wdon1_TVHS_IV_eventplot_RF_FS_robust1.gph, replace
graph export ${figures}/wdon1_TVHS_IV_eventplot_RF_FS_robust1.pdf, as(pdf) replace



*[FIGURE 3] Participation, with decomposition
replace Shh_k3371_dummy5 = 0 if Shh_k3371_dummy_xb==2012
replace Shh_k3371_dummy6 = 0 if Shh_k3371_dummy_xb==2012
replace Shh_k3371_cmlexit5 = 0 if Shh_k3371_cmlexit_xb==2012
replace Shh_k3371_cmlexit6 = 0 if Shh_k3371_cmlexit_xb==2012

#delimit;
twoway  
	(bar Shh_k3371_cmlexit1  	Shh_k3371_cmlexit_xb ,	color(gs3%25)    	  yaxis(1) lpattern(solid) lwidth(medthick) lcolor(none) barwidth(0.5) )
	(bar Shh_k3371_cmlentry1  	Shh_k3371_cmlentry_xb ,	color(blue%25)    	  yaxis(1) lpattern(solid) lwidth(medthick) lcolor(none) barwidth(0.5) )
	(connected Shh_k3371_dummy1 		Shh_k3371_dummy_xb ,	color(blue)     msymbol(s) lpattern(solid) lwidth(medthick))
	(line Shh_k3371_dummy5 		Shh_k3371_dummy_xb ,	color(blue%50)    yaxis(1) msymbol(s) lpattern(longdash) lwidth(medium))
	(line Shh_k3371_dummy6		Shh_k3371_dummy_xb ,	color(blue%50)    yaxis(1) msymbol(s) lpattern(longdash) lwidth(medium))
	, graphregion(color(white))
	ylabel(-0.006(0.002)0.006, axis(1) gmin gmax glcolor(gs14))
	
	xlabel(2006(2)2018)
	xline(2012, lcolor(gs4%15) lwidth(vvthick))
	ytitle("Probability of Giving", axis(1) color(blue))
	legend(order(3 "Effect on whether someone gives" - "  Decomposed into:" 1 "  Exit from giving" 2 "  Entry into giving" ) ring(0) pos(11) 
			region(color(none) lcolor(none)) margin(r=8) bmargin(1 0.5 0.5 5) cols(1) symxsize(*0.5))
	xtitle("Year")
	xsize(1.75) ysize(1);
#delimit cr
graph save ${figures}/wdon1_TVHS_IV_eventplot_extmargin_exit_entry.gph, replace
graph export ${figures}/wdon1_TVHS_IV_eventplot_extmargin_exit_entry.pdf, as(pdf) replace


*[FIGURE A.18] log(Giving) separately for ex-ante givers and nongivers
replace Shh_logk3371_pregiv5 = 0 if Shh_logk3371_pregiv_xb==2012
replace Shh_logk3371_pregiv6 = 0 if Shh_logk3371_pregiv_xb==2012
replace Shh_logk3371_prenongiv5 = 0 if Shh_logk3371_prenongiv_xb==2012
replace Shh_logk3371_prenongiv6 = 0 if Shh_logk3371_prenongiv_xb==2012
#delimit;
twoway  
	(connected Shh_logk33711 	Shh_logk3371_xb 	,		lcolor(gs2%25)   mlcolor(none) yaxis(1) msymbol(o) mcolor(gs12) lwidth(medium) lpattern(solid))
	(connected Shh_logk3371_pregiv1 Shh_logk3371_pregiv_xb ,		lcolor(blue)    yaxis(1) msymbol(s) mcolor(blue%75) lwidth(medthick))
	(line Shh_logk3371_pregiv6  	Shh_logk3371_pregiv_xb ,		color(blue%50) yaxis(1) lpattern(longdash))
	(line Shh_logk3371_pregiv5  	Shh_logk3371_pregiv_xb ,		color(blue%50) yaxis(1) lpattern(longdash))
	, graphregion(color(white))
	ylabel(-0.02(0.005)0.005, axis(1) gmin gmax glcolor(gs14))
	
	xlabel(2006(2)2018)
	xmlabel(2007(1)2018, nolabel)
	xline(2012, lcolor(gs4%15) lwidth(vthick))
	ytitle("log of Charitable Giving (2012=0)", axis(1) color(blue))
	
	legend(order(2 "Conditional on giving in pre-period" 1 "Full sample (baseline)" ) ring(0) pos(7) 
			region(color(none) lcolor(none)) margin(r=8) bmargin(1 0.5 1.5 5) cols(1) symxsize(*1))
	xtitle("Year")
	xsize(1.75) ysize(1);
#delimit cr
graph save ${figures}/wdon1_TVHS_IV_eventplot_logk3371_pregiv.gph, replace
graph export ${figures}/wdon1_TVHS_IV_eventplot_logk3371_pregiv.pdf, as(pdf) replace


*[FIGURE A.14] Total giving w/w/o winsorizing
gen Shh_gaver_tot_w_xb_shifted = Shh_gaver_tot_w_xb - 0.15
replace Shh_gaver_tot_w1 	= . if missing(Shh_gaver_tot_w_xb)
replace Shh_gaver_tot_w5 	= . if missing(Shh_gaver_tot_w_xb)
replace Shh_gaver_tot_w6 	= . if missing(Shh_gaver_tot_w_xb)
#delimit;
twoway  
	(connected Shh_gaver_tot1 		Shh_gaver_tot_xb,		color(blue)    yaxis(1) msymbol(Sh))
	(rcap Shh_gaver_tot6 Shh_gaver_tot5 	Shh_gaver_tot_xb,		color(blue%50) yaxis(1))
	(scatter Shh_gaver_tot_w1 		Shh_gaver_tot_w_xb_shifted ,	color(red)    yaxis(1) msymbol(Oh) lpattern(shortdash))
	(rcap Shh_gaver_tot_w6 Shh_gaver_tot_w5 Shh_gaver_tot_w_xb_shifted ,	color(red%50) yaxis(1))
	
	, graphregion(color(white))
	ylabel(-60(20)60, axis(1) gmin gmax glcolor(gs14))
	
	xlabel(2011(2)2018)
	xline(2012, lcolor(gs4%15) lwidth(vvvthick))
	ytitle("Amount of Charitable Giving (NOK, 2012=0)", axis(1))
	legend(order(1 "Giving (raw)" 3  "Giving (winsorized)"  ) ring(0) pos(11) region(lcolor(gs12)) margin(r=8) bmargin(1 0.5 0.5 4) cols(1))
	xtitle("Year")
	xsize(1.75) ysize(1);
#delimit cr
graph save ${figures}/wdon1_robust_winsorizing_1.gph, replace
graph export ${figures}/wdon1_robust_winsorizing_1.pdf, as(pdf) replace

*[FIGURE A.16] log(Giving) effect by frequency of giving during 2006-2012
replace Shh_logk3371_regular5 = 0 if Shh_logk3371_regular_xb==2012
replace Shh_logk3371_regular6 = 0 if Shh_logk3371_regular_xb==2012
replace Shh_logk3371_occasional5 = 0 if Shh_logk3371_occasional_xb==2012
replace Shh_logk3371_occasional6 = 0 if Shh_logk3371_occasional_xb==2012
#delimit;
twoway  
	(connected 	Shh_logk3371_occasional1 	Shh_logk3371_occasional_xb 	,	lcolor(green)   mlcolor(none) yaxis(1) msymbol(O) mcolor(green) lwidth(medium) lpattern(solid))
	(line 		Shh_logk3371_occasional6  	Shh_logk3371_occasional_xb ,		color(green%50) yaxis(1) lpattern(shortdash))
	(line 		Shh_logk3371_occasional5  	Shh_logk3371_occasional_xb ,		color(green%50) yaxis(1) lpattern(shortdash))
	(connected 	Shh_logk3371_regular1 		Shh_logk3371_regular_xb ,		lcolor(blue)    yaxis(1) msymbol(s) mcolor(blue%75) lwidth(medthick))
	(line 		Shh_logk3371_regular6  		Shh_logk3371_regular_xb ,		color(blue%50) yaxis(1) lpattern(longdash))
	(line 		Shh_logk3371_regular5  		Shh_logk3371_regular_xb ,		color(blue%50) yaxis(1) lpattern(longdash))
	, graphregion(color(white))
	ylabel(-0.025(0.005)0.005, axis(1) gmin gmax glcolor(gs14))
	
	xlabel(2012(2)2018)
	xmlabel(2012(1)2018, nolabel)
	xline(2012, lcolor(gs4%15) lwidth(vthick))
	ytitle("log of Charitable Giving (2012=0)", axis(1) color(black))
	
	legend(order(4 "Gave regularly (6 or 7 out of 7 pre-period years)"  1 "Gave occasionally (1, 2, or 3 per-period years)"  ) ring(0) pos(7) 
			region(color(none) lcolor(none)) margin(r=8) bmargin(1 0.5 1.5 5) cols(1) symxsize(*1))
	xtitle("Year")
	xsize(1.33) ysize(1);
#delimit cr
graph save ${figures}/wdon1_robust_event_logk3371_occasional_regular.gph, replace
graph export ${figures}/wdon1_robust_event_logk3371_occasional_regular.pdf, as(pdf) replace


*-----------------------------------------------*
* IV: Instrument interacted with TNW bins
*-----------------------------------------------*
	preserve
	keep if $if_TVHS_cTNW & hh_head_2012 ==1 & inrange(year,2012,2018)
	foreach Yvar in    "log_g_regular" "log_g_occasional"  "log_g_liquid" "log_g" "log_g_dom_c12"  "log_g_rel_c12" "log_g_nonrel_c12" "log_g_intl_c12" "log_g_condpos" "g_pos" "g_pos_prenongiv"  "log_g_pregiv"  "log_g" "hh_gaver_tot_w"  "FD_log_g" {
		global Yvar = "`Yvar'"
		global slist = "1 2 3"
		global regname = "$Yvar"
		if "$Yvar" == "hh_gaver_tot_w" {
			global regname = "amount"
			global slist = "1 2 3"
		}
		if "$Yvar" == "FD_log_g" {
			global regname = "FD_log_g"
			global slist = "4 5 6"
		}
		if "$Yvar" == "log_g" {
			global regname = "log_g"
			global slist = "4 5 6"
		}
		if inlist("$Yvar","log_g_pregiv","g_pos_prenongiv","g_pos","log_g_condpos","log_g_rel","log_g_nonrel","log_g_intl","log_g_dom") {
			global slist = "4"
		}
		if inlist("$Yvar","log_g_rel_c12", "log_g_nonrel_c12", "log_g_intl_c12","log_g_dom_c12") {
			global slist = "4"
		}
		if inlist("$Yvar","log_g_regular", "log_g_occasional") {
			global slist = "4"
			global regname = subinstr("$Yvar","log_","l",.)
		}
		if "$Yvar" == "log_g_liquid" {
			global regname = "log_g_liquid"
			global slist = "4"
		}
		foreach s of numlist  $slist {
			global s = `s'
			#delimit;
			if $s == 1 {; local xvarlist = "hh_wtax_dummy hh_wtax_b12"; 	};
			if $s == 2 {; local xvarlist = "hh_wtax_b12"; 			};
			if $s == 3 {; local xvarlist = "hh_wtax_dummy"; 		};
			if $s == 4 {; local xvarlist = "hh_AWTR hh_MWTR"; 		};
			if $s == 5 {; local xvarlist = "hh_AWTR"; 			};
			if $s == 6 {; local xvarlist = "hh_MWTR"; 			};
			if $s == 7 {; local xvarlist = "D.hh_AWTR D.hh_MWTR"; 		};
			#delimit cr
			global xvarlist = "`xvarlist'"
			
			disp "[ $S_DATE $S_TIME ] <------------> start IV  <------------>  Yvar = ${Yvar}, s=${s},  slist =${slist}, xvarlist = ${xvarlist}, regname = ${regname}"
			#delimit;
				ivreghdfe $Yvar
					(  $xvarlist = 
						_hh_ETNW_2012_bin_100k#c.hh_MVH_secM_2012#c.not2012
					)
					if $if_TVHS_cTNW_2 & hh_head_2012==1 & inrange(year,2012,2018)
					& _hh_ETNW_2012_bin_100k >=995
					, cluster(hhid)
					absorb(year##c.(	hh_TNW_2012	hh_TNW_2012_sq 	hh_TNW_2012_cu
								age age2 age3 famsize famsize2  hh_loggrossinc_2012 
								hh_N_wtaxpayers
								hh_MVH_cabsec_dummy_2012
								
							)
						
						(year#_hh_ETNW_2012_bin_100k)##c.(hh_MVH_priseccabM_2012) 
						_hh_ETNW_2012_bin_100k##c.hh_MVH_secM_2012
						)
					saverf savefirst rf first
					;
				#delimit cr
			
			
			if "$xvarlist" == "hh_AWTR hh_MWTR" {
				nlcom _b[hh_AWTR] + _b[hh_MWTR]
				estadd scalar SUM_b = _b[hh_AWTR] + _b[hh_MWTR]
				estadd scalar SUM_se = sqrt(r(V)[1,1])
			}

			estimates save ${temp_estimates}/wdon1_${regname}_TVHS_IV_${s}_bybins, replace
			estimates restore _ivreg2_${Yvar}
			estimates save ${temp_estimates}/wdon1_${regname}_TVHS_IV_${s}_RF_bybins, replace
			
			* store first-stage regs. extensive-marg stuff (wtaxdum, mwtr) is FS1; intensive is FS2
			foreach xvar in $xvarlist {
				estimates restore _ivreg2_`xvar'
				if "`xvar'"=="hh_wtax_dummy" | "`xvar'"=="hh_MWTR" {
					estimates save ${temp_estimates}/wdon1_${regname}_TVHS_IV_${s}_FS1_bybins, replace
				}
				if "`xvar'"=="hh_wtax_b12" | "`xvar'"=="hh_AWTR"{
					estimates save ${temp_estimates}/wdon1_${regname}_TVHS_IV_${s}_FS2_bybins, replace
				}
				
			}
		}
	}
	restore
	

*------------------------*
* Prep tables
*------------------------*
*--------------------------------------------------*		
* IV Table : instrument # TNW bins
*--------------------------------------------------*	
* Load all IV regressions for use with esttab
global regname = "amount"
foreach X in 	"wdon1_${regname}_TVHS_IV_1_bybins" "wdon1_${regname}_TVHS_IV_1_RF_bybins"	"wdon1_${regname}_TVHS_IV_1_FS1_bybins"	"wdon1_${regname}_TVHS_IV_1_FS2_bybins"	///
		"wdon1_${regname}_TVHS_IV_2_bybins" "wdon1_${regname}_TVHS_IV_2_RF_bybins"				  		"wdon1_${regname}_TVHS_IV_1_FS2_bybins"	///
		"wdon1_${regname}_TVHS_IV_3_bybins" "wdon1_${regname}_TVHS_IV_3_RF_bybins"	"wdon1_${regname}_TVHS_IV_1_FS1_bybins"	{
	estimates use ${temp_estimates}/`X'
	local eststo_name = subinstr("`X'","wdon1_","",.) /* shorten name for eststo command compatibility*/
	local eststo_name = subinstr("`eststo_name'","_TVHS","",.)
	eststo `eststo_name': estimates replay	
}
global regname = "log_g"
foreach X in 	"wdon1_${regname}_TVHS_IV_4_bybins" "wdon1_${regname}_TVHS_IV_4_RF_bybins"	"wdon1_${regname}_TVHS_IV_4_FS1_bybins"	"wdon1_${regname}_TVHS_IV_4_FS2_bybins"	///
		"wdon1_${regname}_TVHS_IV_5_bybins" "wdon1_${regname}_TVHS_IV_5_RF_bybins"				  		"wdon1_${regname}_TVHS_IV_5_FS2_bybins"	///
		"wdon1_${regname}_TVHS_IV_6_bybins" "wdon1_${regname}_TVHS_IV_6_RF_bybins"	"wdon1_${regname}_TVHS_IV_6_FS1_bybins"	 					 {
	estimates use ${temp_estimates}/`X'
	local eststo_name = subinstr("`X'","wdon1_","",.)
	local eststo_name = subinstr("`eststo_name'","_TVHS","",.)
	eststo `eststo_name': estimates replay	
}
global regname = "FD_log_g"
foreach X in 	"wdon1_${regname}_TVHS_IV_4_bybins" "wdon1_${regname}_TVHS_IV_4_RF_bybins"	"wdon1_${regname}_TVHS_IV_4_FS1_bybins"	"wdon1_${regname}_TVHS_IV_4_FS2_bybins"	///
		"wdon1_${regname}_TVHS_IV_5_bybins" "wdon1_${regname}_TVHS_IV_5_RF_bybins"				  		"wdon1_${regname}_TVHS_IV_5_FS2_bybins"	///
		"wdon1_${regname}_TVHS_IV_6_bybins" "wdon1_${regname}_TVHS_IV_6_RF_bybins"	"wdon1_${regname}_TVHS_IV_6_FS1_bybins"	{
	estimates use ${temp_estimates}/`X'
	local eststo_name = subinstr("`X'","wdon1_","",.)
	local eststo_name = subinstr("`eststo_name'","_TVHS","",.)
	eststo `eststo_name': estimates replay	
}

global regname = "log_g_condpos"
foreach X in 	"wdon1_${regname}_TVHS_IV_4_bybins" "wdon1_${regname}_TVHS_IV_4_RF_bybins"	"wdon1_${regname}_TVHS_IV_4_FS1_bybins"	"wdon1_${regname}_TVHS_IV_4_FS2_bybins"	{
	estimates use ${temp_estimates}/`X'
	local eststo_name = subinstr("`X'","wdon1_","",.)
	local eststo_name = subinstr("`eststo_name'","condpos","cp",.)
	local eststo_name = subinstr("`eststo_name'","_TVHS","",.)
	eststo `eststo_name': estimates replay	
}

global regname = "g_pos_prenongiv"
foreach X in 	"wdon1_${regname}_TVHS_IV_4_bybins" "wdon1_${regname}_TVHS_IV_4_RF_bybins"	"wdon1_${regname}_TVHS_IV_4_FS1_bybins"	"wdon1_${regname}_TVHS_IV_4_FS2_bybins"	{
	estimates use ${temp_estimates}/`X'
	local eststo_name = subinstr("`X'","wdon1_","",.)
	local eststo_name = subinstr("`eststo_name'","prenongiv","png",.)
	local eststo_name = subinstr("`eststo_name'","_TVHS","",.)
	eststo `eststo_name': estimates replay	
}


global regname = "log_g_liquid"
foreach X in 	"wdon1_${regname}_TVHS_IV_4_bybins" "wdon1_${regname}_TVHS_IV_4_RF_bybins"	"wdon1_${regname}_TVHS_IV_4_FS1_bybins"	"wdon1_${regname}_TVHS_IV_4_FS2_bybins"	{
	estimates use ${temp_estimates}/`X'
	local eststo_name = subinstr("`X'","wdon1_","",.)
	local eststo_name = subinstr("`eststo_name'","prenongiv","png",.)
	local eststo_name = subinstr("`eststo_name'","_TVHS","",.)
	eststo `eststo_name': estimates replay	
}

foreach var of varlist log_g_rel_c12 log_g_nonrel_c12 log_g_intl_c12 log_g_dom_c12 {
	global regname = subinstr("`var'","log_","l",.)
	foreach X in "wdon1_${regname}_TVHS_IV_4_bybins" "wdon1_${regname}_TVHS_IV_4_RF_bybins"	"wdon1_${regname}_TVHS_IV_4_FS1_bybins"	"wdon1_${regname}_TVHS_IV_4_FS2_bybins"	{
		estimates use ${temp_estimates}/`X'
		local eststo_name = subinstr("`X'","wdon1_","",.)
		local eststo_name = subinstr("`eststo_name'","_TVHS","",.)
		 eststo `eststo_name': estimates replay	
		
		* add some sumstats
		local gvar = subinstr(subinstr("`var'","log_","",.),"_c12","",.)
		capture drop _dum_pos
		gen _dum_pos = `gvar'>0 if !missing(`gvar')
		sum _dum_pos if $if_TVHS_cTNW & hh_head_2012==1 & inrange(year,2012,2018) & hh_MVH_prisecM_2012>0
		estadd scalar ProbPos = `r(mean)' : `eststo_name'
		capture drop _wins
		gen _wins = min(`gvar',100000) if !missing(`gvar')
		sum _wins if $if_TVHS_cTNW & hh_head_2012==1 & inrange(year,2012,2018) & hh_MVH_prisecM_2012>0 & _wins>0 
		estadd scalar CondMean = `r(mean)' : `eststo_name'
	}	
}

*[Table 1] Wealth taxation and charitable giving: main results...
#delimit; 
esttab  amount_IV_2_bybins amount_IV_1_bybins
	log_g_IV_6_bybins log_g_IV_4_bybins
	log_g_cp_IV_4_bybins
	g_pos_png_IV_4_bybins
	, 
	b(5) se(5)
	varwidth(30) stats(SUM_b SUM_se rkf r2 N, fmt(3 3 2 4 0 )) 
	star(* 0.05 ** 0.01 *** 0.001)  nocons
	substitute ( 	#c.not2012 "$ \times \mathds{1}[t>2012]$"
			
			hh_wtax_dummy "1[wtax$ _{h,t}>$0]"
			hh_wtax_b12 "wtax$ _{h,t}$ (NOK)"
			hh_AWTR "AWTR"
			hh_MWTR "MWTR"
			rkf "rk-$ F$-statistic"
			_hh_ETNW_2012_bin_100k ""
			c. ""
			co. ""
			#  "$ \times$"
			hh_MVH_secM_2012 "$ MVHS_{h,12}$"
			
			985b. "1[b=-1.5M]]"
			986.  "1[b=-1.4M]"
			987.  "1[b=-1.3M]"
			988.  "1[b=-1.2M]"
			989.  "1[b=-1.1M]"
			990.  "1[b=-1.0M]"
			991.  "1[b=-0.9M]"
			992.  "1[b=-0.8M]"
			993.  "1[b=-0.7M]"
			994.  "1[b=-0.6M]"
			995b.  "1[b=-0.5M]"
			995.  "1[b=-0.5M]"
			996.  "1[b=-0.4M]"
			997.  "1[b=-0.3M]"
			998.  "1[b=-0.2M]"
			999.  "1[b=-0.1M]"
			1000. "1[b= 0.0M]"
			1001. "1[b= 0.1M]"
			1002. "1[b= 0.2M]"
			1003. "1[b= 0.3M]"
			1004. "1[b= 0.4M]"
			1006. "1[b= 0.6M]"
			1008. "1[b= 0.8M]"
			1012. "1[b= 1.2M]"
			1log_g_r016. "1[b= 1.6M]"
			1020. "1[b= 2.0M]"
			
			
	)
	booktabs nolines nonum nomtitle  frag replace;
#delimit cr


*[TABLE A.8] Robustness: Effect of wealth taxation on charitable giving when removing illiquid households
#delimit; 
esttab  log_g_liquid_IV_4_bybins
	, 
	b(5) se(5)
	varwidth(30) stats(SUM_b SUM_se rkf r2 N, fmt(3 3 2 4 0 )) 
	star(* 0.05 ** 0.01 *** 0.001)  nocons
	substitute ( 	#c.not2012 "$ \times \mathds{1}[t>2012]$"
			
			hh_wtax_dummy "1[wtax$ _{h,t}>$0]"
			hh_wtax_b12 "wtax$ _{h,t}$ (NOK)"
			hh_AWTR "AWTR"
			hh_MWTR "MWTR"
			rkf "rk-$ F$-statistic"
			_hh_ETNW_2012_bin_100k ""
			c. ""
			co. ""
			#  "$ \times$"
			hh_MVH_secM_2012 "$ MVHS_{h,12}$"
			
			
			
	)
	booktabs nolines nonum nomtitle  frag replace;
#delimit cr

*[TABLE A.4] OUTPUT Main Table's underlying REDUCED-FORM coefficients
#delimit; 
esttab  
	amount_IV_1_RF_bybins 
	log_g_IV_4_RF_bybins 
	log_g_cp_IV_4_RF_bybins
	g_pos_png_IV_4_RF_bybins
, 
	b(4) se(4)
	varwidth(30) stats(SUM_b SUM_se rkf r2 N, fmt(3 3 2 4 0 )) 
	/*star(* 0.05 ** 0.01 *** 0.001)*/  nocons
	star(* 0.10 ** 0.05 *** 0.01)
	substitute ( 	#c.not2012 "$ \times \mathds{1}[t>2012]$"
			
			hh_wtax_dummy "1[wtax$ _{h,t}>$0]"
			hh_wtax_b12 "wtax$ _{h,t}$ (NOK)"
			hh_AWTR "AWTR"
			hh_MWTR "MWTR"
			rkf "rk-$ F$-statistic"
			_hh_ETNW_2012_bin_100k ""
			c. ""
			co. ""
			#  "$ \times$"
			hh_MVH_secM_2012 "$ MVHS_{h,12}$"
			
			985b. "1[b=-1.5M]]"
			986.  "1[b=-1.4M]"
			987.  "1[b=-1.3M]"
			988.  "1[b=-1.2M]"
			989.  "1[b=-1.1M]"
			990.  "1[b=-1.0M]"
			991.  "1[b=-0.9M]"
			992.  "1[b=-0.8M]"
			993.  "1[b=-0.7M]"
			994.  "1[b=-0.6M]"
			995b.  "1[b=-0.5M]"
			995.  "1[b=-0.5M]"
			996.  "1[b=-0.4M]"
			997.  "1[b=-0.3M]"
			998.  "1[b=-0.2M]"
			999.  "1[b=-0.1M]"
			1000. "1[b= 0.0M]"
			1001. "1[b= 0.1M]"
			1002. "1[b= 0.2M]"
			1003. "1[b= 0.3M]"
			1004. "1[b= 0.4M]"
			1006. "1[b= 0.6M]"
			1008. "1[b= 0.8M]"
			1012. "1[b= 1.2M]"
			1016. "1[b= 1.6M]"
			1020. "1[b= 2.0M]"
			
			
	)
	booktabs nolines nonum nomtitle  frag replace;
#delimit cr


*[TABLE A.5] OUTPUT Main Table's underlying FIRST-STAGE coefficients
#delimit; 
esttab  amount_IV_1_FS1_bybins amount_IV_1_FS2_bybins
	log_g_IV_4_FS1_bybins log_g_IV_4_FS2_bybins
	
,
	b(5) se(5)
	varwidth(30) stats(SUM_b SUM_se rkf r2 N, fmt(3 3 2 4 0 )) 
	star(* 0.05 ** 0.01 *** 0.001)  nocons
	substitute ( 	#c.not2012 "$ \times \mathds{1}[t>2012]$"
			
			hh_wtax_dummy "1[wtax$ _{h,t}>$0]"
			hh_wtax_b12 "wtax$ _{h,t}$ (NOK)"
			hh_AWTR "AWTR"
			hh_MWTR "MWTR"
			rkf "rk-$ F$-statistic"
			_hh_ETNW_2012_bin_100k ""
			c. ""
			co. ""
			#  "$ \times$"
			hh_MVH_secM_2012 "$ MVHS_{h,12}$"
			
			985b. "1[b=-1.5M]]"
			986.  "1[b=-1.4M]"
			987.  "1[b=-1.3M]"
			988.  "1[b=-1.2M]"
			989.  "1[b=-1.1M]"
			990.  "1[b=-1.0M]"
			991.  "1[b=-0.9M]"
			992.  "1[b=-0.8M]"
			993.  "1[b=-0.7M]"
			994.  "1[b=-0.6M]"
			995b.  "1[b=-0.5M]"
			995.  "1[b=-0.5M]"
			996.  "1[b=-0.4M]"
			997.  "1[b=-0.3M]"
			998.  "1[b=-0.2M]"
			999.  "1[b=-0.1M]"
			1000. "1[b= 0.0M]"
			1001. "1[b= 0.1M]"
			1002. "1[b= 0.2M]"
			1003. "1[b= 0.3M]"
			1004. "1[b= 0.4M]"
			1006. "1[b= 0.6M]"
			1008. "1[b= 0.8M]"
			1012. "1[b= 1.2M]"
			1016. "1[b= 1.6M]"
			1020. "1[b= 2.0M]"
			
			
	)
	booktabs nolines nonum nomtitle  frag replace;
#delimit cr


*[TABLE A.7] Table: Religious, Nonrel, Intl
#delimit; 
esttab  lg_rel_c12_IV_4_bybins
	lg_nonrel_c12_IV_4_bybins
	lg_intl_c12_IV_4_bybins
	lg_dom_c12_IV_4_bybins
, 
	b(2) se(2)
	varwidth(30) stats(SUM_b SUM_se ProbPos CondMean rkf N, fmt(2 2 2 0 2 0 )) 
	star(* 0.05 ** 0.01 *** 0.001)  nocons
	substitute ( 	#c.not2012 "$ \times \mathds{1}[t>2012]$"
			
			hh_wtax_dummy "1[wtax$ _{h,t}>$0]"
			hh_wtax_b12 "wtax$ _{h,t}$ (NOK)"
			hh_AWTR "AWTR"
			hh_MWTR "MWTR"
			rkf "rk-$ F$-statistic"
			_hh_ETNW_2012_bin_100k ""
			c. ""
			co. ""
			#  "$ \times$"
			
			
			
	)
	booktabs nolines nonum nomtitle  frag replace;
#delimit cr



*[FIGURE A.1] Figure: Religious, Nonreligious, Intl, .. Giving
local ncounter = -2
foreach var of varlist log_g_rel_c12 log_g_nonrel_c12 log_g_intl_c12 log_g_dom_c12 {
	global regname = subinstr("`var'","log_","l",.)
	local X = "wdon1_${regname}_TVHS_IV_4_bybins" 	
		estimates use ${temp_estimates}/`X'
		local eststo_name = subinstr("`X'","wdon1_","",.)
		local eststo_name = subinstr("`eststo_name'","_TVHS","",.)
		 eststo `eststo_name': estimates replay	
	
	local ncounter = `ncounter'+3
	local ncounter_plus2 = `ncounter'+5
	capture drop _${regname}_*
	gen _${regname}_b = e(SUM_b) 		     in `ncounter'/`ncounter_plus2'
	gen _${regname}_LB = e(SUM_b)-1.96*e(SUM_se) in `ncounter'/`ncounter_plus2'
	gen _${regname}_UB = e(SUM_b)+1.96*e(SUM_se) in `ncounter'/`ncounter_plus2'
}
capture drop _xx 
gen byte _xx = _n in 1/15
global fcol = "green%20"
global lcol = "green"
#delimit;
twoway 	(rarea	_lg_rel_c12_UB _lg_rel_c12_LB _xx		in 1/3,   lcolor(blue) lwidth(none) fcolor(blue%20)	)
	(line	_lg_rel_c12_b 			_xx		in 1/3,   lcolor(blue) lwidth(thick)	lpattern(solid)	)
	(rarea	_lg_nonrel_c12_UB _lg_nonrel_c12_LB	_xx     in 4/6,   lcolor(blue) lwidth(none) fcolor(${fcol})	)
	(line	_lg_nonrel_c12_b 			_xx	in 4/6,   lcolor(${lcol}) lwidth(thick)	lpattern(solid)	)
	(rarea	_lg_intl_c12_UB _lg_intl_c12_LB _xx		in 7/9,   lcolor(blue) lwidth(none) fcolor(green%15)	)
	(line	_lg_intl_c12_b 				_xx	in 7/9,   lcolor(${lcol}) lwidth(medthick)	lpattern(solid)	)
	(rarea	_lg_dom_c12_UB _lg_dom_c12_LB	_xx     	in 10/12, lcolor(blue) lwidth(none) fcolor(green%15)	)
	(line	_lg_dom_c12_b 			_xx		in 10/12, lcolor(${lcol}) lwidth(medthick)	lpattern(solid) )
	, graphregion(color(white)) plotregion(lcolor(none))
	 yscale(noline) xscale(noline) ylabel(,noticks) xlabel(,noticks)
	xline(6.5, lcolor(gs12) lwidth(medthick) lpattern(longdash))
	xlabel(  2 "{stSerif:Religious}"  5 "{stSerif:Nonreligious}" 8 "{stSerif:International}" 11 "{stSerif:Domestic}"  , labsize(medlarge)  )
	ylabel(-100(25)50, glcolor(gs14) gmin gmax)
	legend(off)
	ytitle("Semi-elasticity of Giving wrt. wealth tax rate")
	xtitle("")
	;
#delimit cr
 
graph save ${figures}/wdon1_TVHS_IV_wtax_effect_by_type.gph, replace
graph export ${figures}/wdon1_TVHS_IV_wtax_effect_by_type.pdf, as(pdf) replace

*------------------------------------------
* First-stage Figure
*------------------------------------------
estimates use ${temp_estimates}/wdon1_amount_TVHS_IV_1_FS1_bybins
estimates replay
mat FS_wtax_dummy = r(table)'
estimates use ${temp_estimates}/wdon1_amount_TVHS_IV_1_FS2_bybins
estimates replay
mat FS_wtax = r(table)'
estimates use ${temp_estimates}/wdon1_log_g_TVHS_IV_4_FS1_bybins
estimates replay
mat FS_MWTR = r(table)'
estimates use ${temp_estimates}/wdon1_log_g_TVHS_IV_4_FS2_bybins
estimates replay
mat FS_AWTR = r(table)'


foreach X in  "wtax" "MWTR" "wtax_dummy" "AWTR" {
		global X = "`X'"
		capture drop FS_${X}*
		svmat2  FS_${X}, rnames(FS_${X}_rn) 
		replace FS_${X}_rn = substr(FS_${X}_rn,1,4)
		replace FS_${X}_rn = subinstr(FS_${X}_rn,"b","",.)
		replace FS_${X}_rn = subinstr(FS_${X}_rn,".","",.)
		
		gen FS_${X}_x = (real(FS_${X}_rn)-1000)
		replace FS_${X}_x = FS_${X}_x/10 /* to get in MNOKs */
		drop FS_${X}2 FS_${X}3 FS_${X}4  FS_${X}7 FS_${X}8 FS_${X}9
}
replace FS_MWTR1 = FS_MWTR1 * 100
replace FS_AWTR1 = FS_AWTR1 * 100
replace FS_wtax1 = FS_wtax1/1000


*[FIGURE A.6] First-stage heterogeneity
#delimit;
twoway  
	(connected FS_MWTR1 FS_MWTR_x,			yaxis(1) lcolor(green) mcolor(green)    mlcolor(none) msymbol(o) lwidth(medium) lpattern(dash))
	(connected FS_AWTR1 FS_AWTR_x,			yaxis(1) lcolor(blue) mcolor(blue)    mlcolor(none) msymbol(s) lwidth(medium) lpattern(shortdash))
	(connected FS_wtax_dummy1 FS_wtax_dummy_x,	yaxis(1) lcolor(green) mcolor(green)    mlcolor(none) msymbol(s) lwidth(medium) lpattern(solid))
	(connected FS_wtax1 FS_wtax_x,			yaxis(2) lcolor(blue) mcolor(blue)    mlcolor(none) msymbol(o) lwidth(medium) lpattern(solid))
	, graphregion(color(white))
	ylabel(-0(0.02)0.1, axis(1) gmin gmax glcolor(gs14))
	ylabel(0(1)5, nogrid axis(2))
	
	ytitle("Tax Rates (%), Probability", axis(1) )
	ytitle("Wealth Tax Amount in NOK 1,000s", axis(2))
	legend(order(3 "1[Wealth Tax>0]" 4 "Wealth Tax (NOK 1,000s)" 1 "MWTR (%)"    2 "AWTR (%)"  ) ring(0) pos(2) 
			region(color(none) lcolor(none)) margin( l=1 r=3 t=1) bmargin(1 0.5 1.5 5) rows(2) cols(2) symxsize(*0.75))
	xlabel(-0.5 "-0.5 MNOK" 0 "0" 0.5 "0.5" 1.5 "1.5" 2 "2+ MNOK")
	xtitle("Taxable Net Wealth in Excess of Threshold (2012, MNOK)")
	xsize(1.75) ysize(1);
#delimit cr
graph save ${figures}/wdon1_TVHS_IV_FS_plot.gph, replace
graph export ${figures}/wdon1_TVHS_IV_FS_plot.pdf, as(pdf) replace

* save special version for slides
gr use ${figures}/wdon1_TVHS_IV_FS_plot.gph,
gr_edit .plotregion1.plot1.style.editstyle line(color(blue) width(medthick))  editcopy
gr_edit .plotregion1.plot1.style.editstyle marker(linestyle(color(blue)) fillcolor(blue)	) editcopy
gr_edit .plotregion1.plot3.style.editstyle line(color(blue) width(medthick))  editcopy
gr_edit .plotregion1.plot3.style.editstyle marker(linestyle(color(blue)) fillcolor(blue)	) editcopy
gr_edit .plotregion1.plot2.style.editstyle line(color(red) width(medthick))  editcopy
gr_edit .plotregion1.plot2.style.editstyle marker(linestyle(color(red)) fillcolor(red)	) editcopy
gr_edit .plotregion2.plot4.style.editstyle line(color(red) width(medthick))  editcopy
gr_edit .plotregion2.plot4.style.editstyle marker(linestyle(color(red)) fillcolor(red)	) editcopy
graph export ${figures}/wdon1_TVHS_IV_FS_plot_SLIDES.pdf, as(pdf) replace

 
*===================================================================*
* IV Approach for Elasticity Heterogeneity
*===================================================================*

	preserve
	keep if $ifSRwide & $if_TVHS_cTNW & inrange(year,2012,2018)
	foreach Yvar in  "Ind100" {
		global Yvar = "`Yvar'"
		global regname = "${Yvar}"
		global slist = "5 1 2 3 4 6"
		foreach s of numlist  $slist {
			global s = `s'
			#delimit;
			if $s == 1 {; local xvarlist = "hh_wtax_dummy hh_wtax_b12_NOK1000"; 	};
			if $s == 2 {; local xvarlist = "hh_wtax_b12_NOK1000"; 			};
			if $s == 3 {; local xvarlist = "hh_wtax_dummy"; 		};
			if $s == 4 {; local xvarlist = "hh_AWTR hh_MWTR"; 		};
			if $s == 5 {; local xvarlist = "hh_AWTR"; 			};
			if $s == 6 {; local xvarlist = "hh_MWTR"; 			};
			#delimit cr
			global xvarlist = "`xvarlist'"
			
			disp "[ $S_DATE $S_TIME ] <------------> start IV  <------------>  Yvar = ${Yvar}, s=${s},  slist =${slist}, xvarlist = ${xvarlist}"
			#delimit;
				ivreghdfe $Yvar
						
					(  $xvarlist = 
						_hh_ETNW_2012_bin_wide#c.hh_MVH_secM_2012#c.not2012
					)
					if $ifSRwide & $if_TVHS_cTNW & inrange(year,2012,2018)
					& _hh_ETNW_2012_bin_wide >=996
					, cluster(hhid)
					absorb(year##c.(	
								age age2 age3 famsize famsize2  hh_loggrossinc_2012 
								hh_N_wtaxpayers
								hh_MVH_cabsec_dummy_2012
								_g1 _g2 _g3 /* third-order polynomial in amount of giving */
							)
						
						(year#_hh_ETNW_2012_bin_100k)##c.(hh_MVH_priseccabM_2012)
						_hh_ETNW_2012_bin_100k##c.hh_MVH_secM_2012
						)
						
					saverf savefirst rf first
					;
				#delimit cr
			
			
			if "$xvarlist" == "hh_AWTR hh_MWTR" {
				nlcom _b[hh_AWTR] + _b[hh_MWTR]
				estadd scalar SUM_b = _b[hh_AWTR] + _b[hh_MWTR]
				estadd scalar SUM_se = sqrt(r(V)[1,1])
			}

			estimates save ${temp_estimates}/wdon1_${regname}_TVHS_IV_${s}_bybins, replace
			estimates restore _ivreg2_${Yvar}
			estimates save ${temp_estimates}/wdon1_${regname}_TVHS_IV_${s}_RF_bybins, replace
			
			* store first-stage regs. extensive-marg stuff (wtaxdum, mwtr) is FS1; intensive is FS2
			foreach xvar in $xvarlist {
				estimates restore _ivreg2_`xvar'
				if "`xvar'"=="hh_wtax_dummy" | "`xvar'"=="hh_MWTR" {
					estimates save ${temp_estimates}/wdon1_${regname}_TVHS_IV_${s}_FS1_bybins, replace
				}
				if "`xvar'"=="hh_wtax_b12" | "`xvar'"=="hh_AWTR"{
					estimates save ${temp_estimates}/wdon1_${regname}_TVHS_IV_${s}_FS2_bybins, replace
				}
				
			}
		}
	}
	restore
	
foreach X in 	"wdon1_Ind100_TVHS_IV_2_bybins"  "wdon1_Ind100_TVHS_IV_3_bybins"  "wdon1_Ind100_TVHS_IV_1_bybins"  "wdon1_Ind100_TVHS_IV_4_bybins" {
	estimates use ${temp_estimates}/`X'
	local X = subinstr("`X'","Ind100","Ind",.)
	eststo `X': estimates replay
	
}


*[TABLE A.6] Wealth taxation and price elasticity of giving
#delimit; 
esttab  wdon1_Ind_TVHS_IV_2_bybins  wdon1_Ind_TVHS_IV_3_bybins  wdon1_Ind_TVHS_IV_1_bybins  wdon1_Ind_TVHS_IV_4_bybins
	
	, 
	b(3) se(3)
	varwidth(30) stats(SUM_b SUM_se rkf N, fmt(4 4 2 0 )) 
	 
	star(* 0.05 ** 0.01 *** 0.001)  nocons
	substitute ( 	#c.not2012 ""
			hh_AWTR "AWTR$ _{h,t}$"
			hh_MWTR "MWTR$ _{h,t}$"
			TNWdistcu "$ (T_t-TNW_{i,2012})^3$"
			TNWdistsq "$ (T_t-TNW_{i,2012})^2$"
			TNWdist "$ (T_t-TNW_{i,2012})$"
			hh_loggrossinc "log(Household Gross Income)"
			hh_wtax_dummy "1[wtax$ _{h,t}$ $>$0]"
			hh_wtax_b12 "wtax$ _{h,t}$"
			hh_logMNW "log(Market-value Net Wealth)"
			hh_MVH_secM_2012 "MVHS$ _{h,2012}$"
			c. ""
			_NOK1000 " (NOK 1,000)"
			
			SUM_b "AWTR+MWTR"
			SUM_se ""
			co. ""
			rkf "rk-$ F$-statistic"
			#  "$ \times$"
			wtax_TVHS_rate "$ M_t$"
			hh_MVH_secM_adj_2012 "$ MVHS_{i,2012}$"
	)
	booktabs nolines nonum nomtitle nogaps frag replace;
#delimit cr


*=============================================================================*
*  Robustness Figure: assign treatment in 2010
*=============================================================================*
gen hh_TNW_2010_sq = hh_TNW_2010^2
gen hh_TNW_2010_cu = hh_TNW_2010^3
gen hh_MVH_prisecM_2010 = hh_MVH_priM_2010 + hh_MVH_secM_2010
gen age_2010 = age_2012-2
gen age2_2010 = (age_2012-2)^2
gen age3_2010 = (age_2012-2)^3
gen hh_MVH_cabsec_dummy_2010 = hh_MVH_cabsecM_2010 > 0 if !missing(hh_MVH_cabsecM_2010)


sort LOPNR year
foreach var of varlist famsize hh_loggrossinc hh_N_wtaxpayers hh_MVH_cabsec_dummy hh_head hh_cTNW {
	by LOPNR (year): egen `var'_2010 = total(`var' * (year==2010)), mi
}
gen famsize2_2010 = famsize_2010^2

by LOPNR (year): egen double hhid_2010 = total(hhid * (year==2010)), mi
gen byte not2010 = (year!=2010)
capture drop not2012
gen byte not2012 = (year!=2012)
	global Yvar = "hh_logk3371"
	
	
	
	#delimit;
	reghdfe	$Yvar
			year#c.hh_MVH_secM_2010#c.not2010
			c.hh_MVH_secM_2010
			
		if $if_TVHS_cTNW_2 & hh_head_2012==1
		& inrange(year,2006,2018)
		, cluster(hhid_2010)
		absorb(year##c.(	hh_MVH_priseccabM_2010
					hh_TNW_2010	hh_TNW_2010_sq 	hh_TNW_2010_cu
					age_2010 age2_2010 age3_2010 famsize_2010 famsize2_2010 
					hh_loggrossinc_2010
					hh_N_wtaxpayers_2010
					hh_MVH_cabsec_dummy_2010
				)
			)
		;
	#delimit cr
estimates save ${temp_estimates}/event_${Yvar}_placebo2010, replace
mat T${Yvar} = r(table)'

mat Shh_logk3371 = Thh_logk3371[1..13,1..6]

foreach Yvar of varlist $Yvar {
	
	
	capture drop S`Yvar'*
	svmat2  S`Yvar', rnames(S`Yvar'_rn) 
		
		replace S`Yvar'_rn = substr(S`Yvar'_rn,1,4)
		replace S`Yvar'_rn = subinstr(S`Yvar'_rn,".","",.)
		replace S`Yvar'_rn = subinstr(S`Yvar'_rn,"_","",.)
		replace S`Yvar'_rn = subinstr(S`Yvar'_rn,"b","",.)
		replace S`Yvar'_rn = subinstr(S`Yvar'_rn,"o","",.)
		gen long S`Yvar'_xb = (real(S`Yvar'_rn))
		drop S`Yvar'2 S`Yvar'3  S`Yvar'4 
}

gen Shh_logk3371_xb_shifted = Shh_logk3371_xb - 0.15
replace Shh_logk33715 = Shh_logk33711 if Shh_logk3371_xb==2010
replace Shh_logk33716 = Shh_logk33711 if Shh_logk3371_xb==2010

* for reference: get main estimate
estimates use ${temp_estimates}/event_${Yvar}
estimates replay
mat T${Yvar} = r(table)'

mat Bhh_logk3371 = Thh_logk3371[1..13,1..6]

foreach Yvar of varlist $Yvar {
	
	
	capture drop B`Yvar'*
	svmat2  B`Yvar', rnames(B`Yvar'_rn) 
		
		replace  B`Yvar'_rn =   substr(B`Yvar'_rn,1,4)
		replace  B`Yvar'_rn = subinstr(B`Yvar'_rn,".","",.)
		replace  B`Yvar'_rn = subinstr(B`Yvar'_rn,"_","",.)
		replace  B`Yvar'_rn = subinstr(B`Yvar'_rn,"b","",.)
		replace  B`Yvar'_rn = subinstr(B`Yvar'_rn,"o","",.)
		gen long B`Yvar'_xb =     real(B`Yvar'_rn)
		drop B`Yvar'2 B`Yvar'3  B`Yvar'4 
}

gen Bhh_logk3371_xb_shifted = Bhh_logk3371_xb - 0.15
replace Bhh_logk33715 = Bhh_logk33711 if Bhh_logk3371_xb==2012
replace Bhh_logk33716 = Bhh_logk33711 if Bhh_logk3371_xb==2012

*[FIGURE A.9] Robustness, earlier treatment assignment
#delimit;
twoway  
	(connected 	Shh_logk33711 	Shh_logk3371_xb_shifted if Shh_logk3371_xb >=2005.5,		lcolor(orange)    yaxis(1) msymbol(o) mcolor(orange%75) lwidth(medthick))
	(line 		Shh_logk33716 	Shh_logk3371_xb_shifted if Shh_logk3371_xb >=2005.5,		color(orange%40) yaxis(1) lpattern(shortdash))
	(line		Shh_logk33715 	Shh_logk3371_xb_shifted if Shh_logk3371_xb >=2005.5,		color(orange%40) yaxis(1) lpattern(shortdash))
	
	(connected 	Bhh_logk33711 	Bhh_logk3371_xb 	if Bhh_logk3371_xb >=2005.5,		lcolor(blue)    yaxis(1) msymbol(s) mcolor(blue%75) lwidth(medthick))
	(line 		Bhh_logk33716 	Bhh_logk3371_xb 	if Bhh_logk3371_xb >=2005.5,		color(blue%50) yaxis(1) lpattern(longdash))
	(line		Bhh_logk33715 	Bhh_logk3371_xb 	if Bhh_logk3371_xb >=2005.5,		color(blue%50) yaxis(1) lpattern(longdash))
	
	(connected 	Shh_logk33711 	Shh_logk3371_xb_shifted if Shh_logk3371_xb ==2010,		yaxis(1) msymbol(O) msize(medlarge) mcolor(orange%100) lwidth(medthick))
	(connected 	Bhh_logk33711 	Shh_logk3371_xb		if Shh_logk3371_xb ==2012,		yaxis(1) msymbol(S) msize(medlarge) mcolor(blue%100) lwidth(medthick))
	, graphregion(color(white))
	ylabel(-0.01(0.005)0.005, axis(1) gmin gmax glcolor(gs14))
	
	xlabel(2007(2)2018)
	xmlabel(2007(1)2018, nolabel)
	ytitle("log of Charitable Giving", axis(1) color(black))
	xline(2012, lcolor(gs4%15) lwidth(vvvthick))
	xline(2009.85, lcolor(orange%25) lwidth(vvvthick))
	ytitle("log of Charitable Giving (2012=0)", axis(1) color(blue))
	legend(order(1 "Base year = 2010" 4 "Base year = 2012" ) ring(0) pos(2) region(lcolor(gs12)) margin(r=8) bmargin(1 0.5 0.5 4) cols(1))
	
	xtitle("Year")
	xsize(1.75) ysize(1);
#delimit cr
graph save ${figures}/wdon1_TVHS_IV_eventplot_RF_log_g_placebo2010assignment.gph, replace
graph export ${figures}/wdon1_TVHS_IV_eventplot_RF_log_g_placebo2010assignment.pdf, as(pdf) replace

*=============================================================================*
*  Robustness Figure, number of children
*=============================================================================*

sort hhid_2012 year
by hhid_2012 year: egen hh_numgk = total(antall_barnebarn)
by hhid_2012 year: egen hh_numk = total(antall_barn)
gen hh_numgk_sq = hh_numgk^2
gen hh_numk_sq = hh_numk^2

gen MVH_X_numk = hh_MVH_priseccabM_2012 * hh_numk
gen TNW_X_numk = hh_TNW_2012 * hh_numk

	global Yvar = "hh_logk3371"
	
	
	
	#delimit;
	reghdfe	$Yvar
			year#c.hh_MVH_secM_2012#c.not2012
			c.hh_MVH_secM_2012
			
		if $if_TVHS_cTNW_2 & hh_head_2012==1
		& inrange(year,2006,2018)
		, cluster(hhid_2012)
		absorb(year##c.(	hh_MVH_priseccabM_2012
					hh_TNW_2012	hh_TNW_2012_sq 	hh_TNW_2012_cu
					age_2012 age2_2012 age3_2012 famsize_2012 famsize2_2012 
					hh_loggrossinc_2012
					hh_N_wtaxpayers_2012
					hh_MVH_cabsec_dummy_2012
					hh_numgk hh_numgk_sq
					hh_numk hh_numk_sq
					MVH_X_numk
					TNW_X_numk
				)
			)
		;
	#delimit cr
estimates save ${temp_estimates}/event_${Yvar}_robustness_numk, replace
mat T${Yvar} = r(table)'

mat Shh_logk3371 = Thh_logk3371[1..13,1..6]

foreach Yvar of varlist $Yvar {
	
	
	capture drop S`Yvar'*
	svmat2  S`Yvar', rnames(S`Yvar'_rn) 
		
		replace S`Yvar'_rn = substr(S`Yvar'_rn,1,4)
		replace S`Yvar'_rn = subinstr(S`Yvar'_rn,".","",.)
		replace S`Yvar'_rn = subinstr(S`Yvar'_rn,"_","",.)
		replace S`Yvar'_rn = subinstr(S`Yvar'_rn,"b","",.)
		replace S`Yvar'_rn = subinstr(S`Yvar'_rn,"o","",.)
		gen long S`Yvar'_xb = (real(S`Yvar'_rn))
		drop S`Yvar'2 S`Yvar'3  S`Yvar'4 
}

gen Shh_logk3371_xb_shifted = Shh_logk3371_xb - 0.15
replace Shh_logk33715 = Shh_logk33711 if Shh_logk3371_xb==2012
replace Shh_logk33716 = Shh_logk33711 if Shh_logk3371_xb==2012

* for reference: get main estimate
estimates use ${temp_estimates}/event_${Yvar}
estimates replay
mat T${Yvar} = r(table)'

mat Bhh_logk3371 = Thh_logk3371[1..13,1..6]

foreach Yvar of varlist $Yvar {
	
	
	capture drop B`Yvar'*
	svmat2  B`Yvar', rnames(B`Yvar'_rn) 
		
		replace  B`Yvar'_rn =   substr(B`Yvar'_rn,1,4)
		replace  B`Yvar'_rn = subinstr(B`Yvar'_rn,".","",.)
		replace  B`Yvar'_rn = subinstr(B`Yvar'_rn,"_","",.)
		replace  B`Yvar'_rn = subinstr(B`Yvar'_rn,"b","",.)
		replace  B`Yvar'_rn = subinstr(B`Yvar'_rn,"o","",.)
		gen long B`Yvar'_xb =     real(B`Yvar'_rn)
		drop B`Yvar'2 B`Yvar'3  B`Yvar'4 
}

*[FIGURE A.10]
gen Bhh_logk3371_xb_shifted = Bhh_logk3371_xb - 0.15
replace Bhh_logk33715 = Bhh_logk33711 if Bhh_logk3371_xb==2012
replace Bhh_logk33716 = Bhh_logk33711 if Bhh_logk3371_xb==2012
#delimit;
twoway  
	(connected 	Shh_logk33711 	Shh_logk3371_xb_shifted if Shh_logk3371_xb >=2005.5,		lcolor(orange)    yaxis(1) msymbol(o) mcolor(orange%75) lwidth(medthick))
	(line 		Shh_logk33716 	Shh_logk3371_xb_shifted if Shh_logk3371_xb >=2005.5,		color(orange%40) yaxis(1) lpattern(shortdash))
	(line		Shh_logk33715 	Shh_logk3371_xb_shifted if Shh_logk3371_xb >=2005.5,		color(orange%40) yaxis(1) lpattern(shortdash))
	
	(connected 	Bhh_logk33711 	Bhh_logk3371_xb 	if Bhh_logk3371_xb >=2005.5,		lcolor(blue)    yaxis(1) msymbol(s) mcolor(blue%75) lwidth(medthick))
	(line 		Bhh_logk33716 	Bhh_logk3371_xb 	if Bhh_logk3371_xb >=2005.5,		color(blue%50) yaxis(1) lpattern(longdash))
	(line		Bhh_logk33715 	Bhh_logk3371_xb 	if Bhh_logk3371_xb >=2005.5,		color(blue%50) yaxis(1) lpattern(longdash))

	, graphregion(color(white))
	ylabel(-0.01(0.005)0.005, axis(1) gmin gmax glcolor(gs14))
	
	xlabel(2007(2)2018)
	xmlabel(2007(1)2018, nolabel)
	ytitle("log of Charitable Giving", axis(1) color(black))
	xline(2012, lcolor(gs4%15) lwidth(vvvthick))
	ytitle("log of Charitable Giving (2012=0)", axis(1) color(blue))
	legend(order(1 "w/ Extra controls" 4 "Base specification" ) ring(0) pos(2) region(lcolor(gs12)) margin(r=8) bmargin(1 0.5 0.5 4) cols(1))
	
	xtitle("Year")
	xsize(1.75) ysize(1);
#delimit cr
graph save ${figures}/wdon1_TVHS_IV_eventplot_RF_log_g_robustness_numk.gph, replace
graph export ${figures}/wdon1_TVHS_IV_eventplot_RF_log_g_robustness_numk.pdf, as(pdf) replace

*=============================================================================*
*  Robustness Figure: sample restriction
*=============================================================================*
	global Yvar = "hh_logk3371"
	#delimit;
	reghdfe	$Yvar
			year#c.hh_MVH_secM_2012#c.not2012
			c.hh_MVH_secM_2012
			
		if 
			 abs(hh_cTNW_2012/(750000*hh_N_wtaxpayers_2012)) <= 2
			& hh_MVH_prisecM_2012>0 & !missing(hh_MVH_prisecM_2012) & hh_head_2012==1
		& inrange(year,2006,2018)
		, cluster(hhid_2012)
		absorb(year##c.(	hh_MVH_priseccabM_2012
					hh_TNW_2012	hh_TNW_2012_sq 	hh_TNW_2012_cu
					age_2012 age2_2012 age3_2012 famsize_2012 famsize2_2012 
					hh_loggrossinc_2012
					hh_N_wtaxpayers_2012
					hh_MVH_cabsec_dummy_2012
					hh_numgk hh_numgk_sq
					hh_numk hh_numk_sq
				)
			)
		;
	#delimit cr
estimates save ${temp_estimates}/event_${Yvar}_robustness_sample_v1, replace
mat T${Yvar} = r(table)'
count if e(sample) & year == 2012

mat Shh_logk3371 = Thh_logk3371[1..13,1..6]

foreach Yvar of varlist $Yvar {
	
	
	capture drop S`Yvar'*
	svmat2  S`Yvar', rnames(S`Yvar'_rn) 
		
		replace S`Yvar'_rn = substr(S`Yvar'_rn,1,4)
		replace S`Yvar'_rn = subinstr(S`Yvar'_rn,".","",.)
		replace S`Yvar'_rn = subinstr(S`Yvar'_rn,"_","",.)
		replace S`Yvar'_rn = subinstr(S`Yvar'_rn,"b","",.)
		replace S`Yvar'_rn = subinstr(S`Yvar'_rn,"o","",.)
		gen long S`Yvar'_xb = (real(S`Yvar'_rn))
		drop S`Yvar'2 S`Yvar'3  S`Yvar'4 
}

gen Shh_logk3371_xb_shifted = Shh_logk3371_xb - 0.15
replace Shh_logk33715 = Shh_logk33711 if Shh_logk3371_xb==2012
replace Shh_logk33716 = Shh_logk33711 if Shh_logk3371_xb==2012

* for reference: get main estimate
estimates use ${temp_estimates}/event_${Yvar}
estimates replay
mat T${Yvar} = r(table)'

mat Bhh_logk3371 = Thh_logk3371[1..13,1..6]

foreach Yvar of varlist $Yvar {
	
	
	capture drop B`Yvar'*
	svmat2  B`Yvar', rnames(B`Yvar'_rn) 
		
		replace  B`Yvar'_rn =   substr(B`Yvar'_rn,1,4)
		replace  B`Yvar'_rn = subinstr(B`Yvar'_rn,".","",.)
		replace  B`Yvar'_rn = subinstr(B`Yvar'_rn,"_","",.)
		replace  B`Yvar'_rn = subinstr(B`Yvar'_rn,"b","",.)
		replace  B`Yvar'_rn = subinstr(B`Yvar'_rn,"o","",.)
		gen long B`Yvar'_xb =     real(B`Yvar'_rn)
		drop B`Yvar'2 B`Yvar'3  B`Yvar'4 
}

gen Bhh_logk3371_xb_shifted = Bhh_logk3371_xb - 0.15
replace Bhh_logk33715 = Bhh_logk33711 if Bhh_logk3371_xb==2012
replace Bhh_logk33716 = Bhh_logk33711 if Bhh_logk3371_xb==2012

*[FIGURE A.12] Robustness to changing sample criteria
#delimit;
twoway  
	(connected 	Shh_logk33711 	Shh_logk3371_xb_shifted if Shh_logk3371_xb >=2005.5,		lcolor(orange)    yaxis(1) msymbol(o) mcolor(orange%75) lwidth(medthick))
	(line 		Shh_logk33716 	Shh_logk3371_xb_shifted if Shh_logk3371_xb >=2005.5,		color(orange%40) yaxis(1) lpattern(shortdash))
	(line		Shh_logk33715 	Shh_logk3371_xb_shifted if Shh_logk3371_xb >=2005.5,		color(orange%40) yaxis(1) lpattern(shortdash))
	
	(connected 	Bhh_logk33711 	Bhh_logk3371_xb 	if Bhh_logk3371_xb >=2005.5,		lcolor(blue)    yaxis(1) msymbol(s) mcolor(blue%75) lwidth(medthick))
	(line 		Bhh_logk33716 	Bhh_logk3371_xb 	if Bhh_logk3371_xb >=2005.5,		color(blue%50) yaxis(1) lpattern(longdash))
	(line		Bhh_logk33715 	Bhh_logk3371_xb 	if Bhh_logk3371_xb >=2005.5,		color(blue%50) yaxis(1) lpattern(longdash))

	, graphregion(color(white))
	ylabel(-0.01(0.005)0.005, axis(1) gmin gmax glcolor(gs14))
	
	xlabel(2007(2)2018)
	xmlabel(2007(1)2018, nolabel)
	ytitle("log of Charitable Giving", axis(1) color(black))
	xline(2012, lcolor(gs4%15) lwidth(vvvthick))
	ytitle("log of Charitable Giving (2012=0)", axis(1) color(blue))
	legend(order(1 "Symmetric sample restriction" 4 "Baseline sample" ) ring(0) pos(2) region(lcolor(gs12)) margin(r=8) bmargin(1 0.5 0.5 4) cols(1))
	
	xtitle("Year")
	xsize(1.75) ysize(1);
#delimit cr
graph save ${figures}/wdon1_TVHS_IV_eventplot_RF_log_g_robustness_sample_v1.gph, replace
graph export ${figures}/wdon1_TVHS_IV_eventplot_RF_log_g_robustness_sample_v1.pdf, as(pdf) replace

*=============================================================================*
*  Robustness Figure: log( c + Giving)
*=============================================================================*

foreach c of numlist 100 500 2000 {
	
	estimates use ${temp_estimates}/event_hh_logk3371_c`c'
	estimates replay
	mat T = r(table)'
	mat S_`c' = T[1..13,1..6]
}

local c = 1000
	estimates use ${temp_estimates}/event_hh_logk3371 /*baseline estimate*/
	estimates replay
	mat T = r(table)'
	mat S_`c' = T[1..13,1..6]

foreach c of numlist 100 500 1000 2000 {
	capture drop S_`c'*
	svmat2  S_`c', rnames(S_`c'_rn) 
			
	replace S_`c'_rn = substr(S_`c'_rn,1,4)
	replace S_`c'_rn = subinstr(S_`c'_rn,".","",.)
	replace S_`c'_rn = subinstr(S_`c'_rn,"_","",.)
	replace S_`c'_rn = subinstr(S_`c'_rn,"b","",.)
	replace S_`c'_rn = subinstr(S_`c'_rn,"o","",.)
	gen long S_`c'_xb = (real(S_`c'_rn))
	drop S_`c'2 S_`c'3  S_`c'4 
}


foreach c of numlist 100 500 1000 2000 {
	capture drop S_`c'5_
	gen S_`c'5_ = S_`c'5
	replace S_`c'5_ = 0 if missing(S_`c'5_) & S_`c'_xb==2012
	
	capture drop S_`c'6_
	gen S_`c'6_ = S_`c'6
	replace S_`c'6_ = 0 if missing(S_`c'6_) & S_`c'_xb==2012
}

*[FIGURE A.11] Sensitivity to changing the log-shifting argument
#delimit;
twoway  
	(connected 	S_1001 		S_100_xb  ,		lcolor(orange%50)    yaxis(1) msymbol(t) mcolor(orange%55) lwidth(medthick))
	(line 		S_1006_ 	S_100_xb ,		color(orange%40) yaxis(1) lpattern(shortdash))
	(line		S_1005_ 	S_100_xb ,		color(orange%40) yaxis(1) lpattern(shortdash))
	
	(connected 	S_5001 		S_500_xb  ,		lcolor(orange%70)    yaxis(1) msymbol(d) mcolor(orange%75) lwidth(medthick))
	(line 		S_5006_ 	S_500_xb ,		color(orange%40) yaxis(1) lpattern(shortdash))
	(line		S_5005_ 	S_500_xb ,		color(orange%40) yaxis(1) lpattern(shortdash))
	
	(connected 	S_10001 	S_1000_xb  ,		lcolor(blue%80)    yaxis(1) msymbol(O) mcolor(blue%75) lwidth(medthick))
	(line 		S_10006_ 	S_1000_xb ,		color(blue%40) yaxis(1) lpattern(shortdash))
	(line		S_5005_ 	S_500_xb ,		color(blue%40) yaxis(1) lpattern(shortdash))
	
	(connected 	S_20001 	S_2000_xb  ,		lcolor(orange%90)    yaxis(1) msymbol(s) mcolor(orange%95) lwidth(medthick))
	(line 		S_20006_ 	S_2000_xb ,		color(orange%40) yaxis(1) lpattern(shortdash))
	(line		S_20005_ 	S_2000_xb ,		color(orange%40) yaxis(1) lpattern(shortdash))
	
	, graphregion(color(white))
	ylabel(-0.03(0.01)0.01, axis(1) gmin gmax glcolor(gs14))
	
	xlabel(2007(2)2018)
	xmlabel(2007(1)2018, nolabel)
	ytitle("log of Charitable Giving", axis(1) color(black))
	xline(2012, lcolor(gs4%15) lwidth(vvvthick))
	ytitle("log of Charitable Giving (2012=0)", axis(1) color(black))
	legend(order(10 "c=2000" 7 "c=1000 (baseline)"  4 "c=500"  1 "c=100"  ) ring(0) pos(7) region(lcolor(gs12)) margin(r=8) bmargin(1 0.5 0.5 4) cols(1))
	
	xtitle("Year")
	xsize(1.75) ysize(1);
#delimit cr

graph save ${figures}/wdon1_TVHS_IV_eventplot_RF_log_g_robustness_logshift.gph, replace
graph export ${figures}/wdon1_TVHS_IV_eventplot_RF_log_g_robustness_logshift.pdf, as(pdf) replace

*=============================================================================*
* Event plot by wealth bins
*=============================================================================*

capture drop _bin
gen long _bin = .
replace _bin = -500 	if inrange(hh_cTNW_2012-hh_N_wtaxpayers_2012*750000,	-500000,	-250000	)
replace _bin = -250 	if inrange(hh_cTNW_2012-hh_N_wtaxpayers_2012*750000,	-250000,	0	)
replace _bin = 0	if inrange(hh_cTNW_2012-hh_N_wtaxpayers_2012*750000,	0,		250000	)
replace _bin = 250 	if inrange(hh_cTNW_2012-hh_N_wtaxpayers_2012*750000,	250000,		500000	)
replace _bin = 500 	if inrange(hh_cTNW_2012-hh_N_wtaxpayers_2012*750000,	500000,		1000000	)
replace _bin = 1000 	if inrange(hh_cTNW_2012-hh_N_wtaxpayers_2012*750000,	1000000,	999999999)
replace _bin = 1000+_bin
tab _bin if year==2012 & $if_TVHS_cTNW_2 & hh_head_2012==1

global Yvar = "hh_logk3371"
#delimit;
	reghdfe	$Yvar
			_bin#year#c.hh_MVH_secM_2012#c.not2012
			
			
		if $if_TVHS_cTNW_2 & hh_head_2012==1
		& inrange(year,2006,2018)
		, cluster(hhid_2012)
		absorb(year##c.(	hh_MVH_priseccabM_2012
					hh_TNW_2012	hh_TNW_2012_sq 	hh_TNW_2012_cu
					age_2012 age2_2012 age3_2012 famsize_2012 famsize2_2012 
					hh_loggrossinc_2012
					hh_N_wtaxpayers_2012
					hh_MVH_cabsec_dummy_2012
					
				)
				
				year#_bin##c.hh_MVH_priseccabM_2012
				_bin#c.hh_MVH_secM_2012
			)
		;
	#delimit cr
estimates save ${temp_estimates}/event_${Yvar}_by_ETNW_bins, replace
mat T = r(table)'

capture drop _year
capture drop _beta_*
capture drop _ci*_*
gen int _year = .
local i = 0
foreach t of numlist 2006(1)2018 {
	local i = `i' + 1
	replace _year = 2005 + `i' in `i'
	
	foreach b of numlist -500 -250 0 250 500 1000 {
		local bn = `b' + 1000
		capture gen _beta_`bn'  = .
		replace _beta_`bn' = T[rownumb(T,"`bn'._bin#`t'.year#c.hh_MVH_secM_2012#c.not2012"),1] in `i'
		capture gen _cih_`bn' = .
		replace _cih_`bn' = T[rownumb(T, "`bn'._bin#`t'.year#c.hh_MVH_secM_2012#c.not2012"),6] in `i'
		capture gen _cil_`bn' = .
		replace _cil_`bn' = T[rownumb(T, "`bn'._bin#`t'.year#c.hh_MVH_secM_2012#c.not2012"),5] in `i'
		
		if `t' == 2012 {
			replace _cih_`bn' = 0
			replace _cil_`bn' = 0
		}
		
		
	}
}

gen _year1 = _year-0.25
gen _year2 = _year-0.15
gen _year3 = _year-0.05
gen _year4 = _year+0.05
gen _year5 = _year+0.15
gen _year6 = _year+0.25

*[FIGURE A.13] (Pre-)Treatment Dynamics Within Wealth Bins
#delimit;
twoway  
	(connected 	_beta_500 	_year  ,	lcolor(orange%50) mcolor(orange%55)    	yaxis(1) msymbol(t)  lwidth(medium))
	(line 		_cih_500 	_year ,		color(orange%50) 			yaxis(1) lpattern(shortdash))
	(line		_cil_500 	_year ,		color(orange%50) 			yaxis(1) lpattern(shortdash))
	
	(connected 	_beta_750 	_year  ,	lcolor(orange%60) mcolor(orange%65)    	yaxis(1) msymbol(d)  lwidth(medthick))
	(line 		_cih_750 	_year ,		color(orange%50) 			yaxis(1) lpattern(shortdash))
	(line		_cil_750 	_year ,		color(orange%50) 			yaxis(1) lpattern(shortdash))
	
	(connected 	_beta_1000 	_year  ,	lcolor(orange%70) mcolor(orange%75)    	yaxis(1) msymbol(o)  lwidth(medthick))
	(line 		_cih_1000 	_year ,		color(orange%50) 			yaxis(1) lpattern(shortdash))
	(line		_cil_1000 	_year ,		color(orange%50) 			yaxis(1) lpattern(shortdash))
	
	(connected 	_beta_1250 	_year  ,	lcolor(orange%80) mcolor(orange%85)    	yaxis(1) msymbol(s)  lwidth(medthick))
	(line 		_cih_1250 	_year ,		color(orange%50) 			yaxis(1) lpattern(shortdash))
	(line		_cil_1250 	_year ,		color(orange%50) 			yaxis(1) lpattern(shortdash))
	
	(connected 	_beta_1500 	_year  ,	lcolor(orange%90) mcolor(orange%95)    	yaxis(1) msymbol(T)  lwidth(medthick))
	(line 		_cih_1500 	_year ,		color(orange%50) 			yaxis(1) lpattern(shortdash))
	(line		_cil_1500 	_year ,		color(orange%50) 			yaxis(1) lpattern(shortdash))
	
	(connected 	_beta_2000 	_year  ,	lcolor(orange%100) mcolor(orange%100)    	yaxis(1) msymbol(t)  lwidth(medthick))
	(line 		_cih_2000 	_year ,		color(orange%50) 			yaxis(1) lpattern(shortdash))
	(line		_cil_2000 	_year ,		color(orange%50) 			yaxis(1) lpattern(shortdash))
	, graphregion(color(white))
	ylabel(-0.03(0.01)0.01, axis(1) gmin gmax glcolor(gs14))
	
	xlabel(2007(2)2018)
	xmlabel(2007(1)2018, nolabel)
	ytitle("log of Charitable Giving", axis(1) color(black))
	xline(2012, lcolor(gs4%15) lwidth(vvvthick))
	ytitle("log of Charitable Giving (2012=0)", axis(1) color(black))
	legend(order(1 "[-500,-250)" 4 "[-250,0)" 7 "[0,250)"  10 "[250,500)" 13 "[500,1000)" 16 "1000+" ) ring(0) pos(7) region(lcolor(gs12)) margin(r=8) bmargin(1 0.5 0.5 4) cols(1))
	
	xtitle("Year")
	xsize(1.75) ysize(1);
#delimit cr

graph save ${figures}/wdon1_TVHS_IV_eventplot_by_ETNW_bins.gph, replace
graph export ${figures}/wdon1_TVHS_IV_eventplot_by_ETNW_bins.pdf, as(pdf) replace


*=================================================================================*
* Fraction who sell secondary home per year
*=================================================================================*


gen byte _Sell = 0 if L.hh_TVH_sec>0 & !missing(L.hh_TVH_sec) & !missing(hh_TVH_sec)
replace  _Sell = 1 if L.hh_TVH_sec>0 & !missing(L.hh_TVH_sec) & !missing(hh_TVH_sec) & hh_TVH_sec==0

*[TABLE A.9] Rate of exit of secondary home ownership
tabstat _Sell if $if_TVHS_cTNW_2 & hh_head_2012==1, by(year) statistics(mean)

*=================================================================================*
* Robustness: Propensity Score Weighting and Trimming to Common Support
*=================================================================================*

* Calculate the propensity score
	#delimit;
	probit	hh_MVH_sec_dummy_2012
			 
			hh_MVH_priseccabM_2012
			hh_TNW_2012	hh_TNW_2012_sq 	hh_TNW_2012_cu
			age_2012 age2_2012 age3_2012 famsize_2012 famsize2_2012
			hh_loggrossinc_2012
			hh_N_wtaxpayers_2012
					
					
		if $if_TVHS_cTNW_2 & hh_head_2012==1
		& year==2012
		, cluster(hhid_2012)
		
		;
	#delimit cr
capture drop _prob
predict _prob if e(sample)
capture drop _ps
by LOPNR (year): egen _ps = total(_prob), mi

* Regression weights
/* see Hirano and Imens (2001): Estimation of Causal Effects Using [...] */
capture drop _psweight
gen 	_psweight = 1/_ps 	if hh_MVH_sec_dummy_2012==1 & !missing(_ps)
replace _psweight = 1/(1-_ps) 	if hh_MVH_sec_dummy_2012==0 & !missing(_ps)

*Summary statistics by whether MVHS>0
preserve
keep if  ${if_TVHS_cTNW_2} & !missing(gaver_tot_w) & !missing(hh_wtax_b12) & !missing(hh_MVH_secM_2012) & !missing(age) & !missing(hh_loggrossinc) & year==2012 & hh_head_2012==1
gen hh_MNW_MNOK = hh_MNW/1000000
gen hh_grossinc_MNOK = hh_grossinc/1000000

foreach SUBSET of numlist 0 1 {
	if `SUBSET'==0 {
		global ifstatement = "hh_MVH_secM_2012 ==0"
	}
	if `SUBSET'==1 {
		global ifstatement = "hh_MVH_secM_2012 >0"
	}
	eststo sumstats1_`SUBSET': estpost tabstat hh_gaver_tot_dummy  		[aw=_psweight] if $ifstatement		, statistics(N mean 		)  columns(statistics)
	eststo sumstats2_`SUBSET': estpost tabstat hh_gaver_tot_w 		[aw=_psweight] if $ifstatement 	& gaver_tot_w>0  	, statistics(N mean p25 p50 p75 )  columns(statistics)
	eststo sumstats3_`SUBSET': estpost tabstat hh_wtax_dummy		[aw=_psweight] if $ifstatement		, statistics(N mean 		)  columns(statistics)
	eststo sumstats4_`SUBSET': estpost tabstat hh_wtax_b12			[aw=_psweight] if $ifstatement & hh_wtax_b12>0  	, statistics(N mean p25 p50 p75 )  columns(statistics)
	eststo sumstats5_`SUBSET': estpost tabstat hh_MVH_priseccabM_2012  	[aw=_psweight] if $ifstatement 		, statistics(N mean p25 p50 p75 )  columns(statistics)
	eststo sumstats6_`SUBSET': estpost tabstat hh_MVH_secM_2012 		[aw=_psweight] if $ifstatement		, statistics(N mean p25 p50 p75 )  columns(statistics)
	eststo sumstats7_`SUBSET': estpost tabstat hh_MNW_MNOK 			[aw=_psweight] if $ifstatement		, statistics(N mean p25 p50 p75 )  columns(statistics)
	eststo sumstats8_`SUBSET': estpost tabstat age 				[aw=_psweight] if $ifstatement		, statistics(N mean p25 p50 p75 )  columns(statistics)
	eststo sumstats9_`SUBSET': estpost tabstat hh_grossinc_MNOK 		[aw=_psweight] if $ifstatement		, statistics(N mean p25 p50 p75 )  columns(statistics)
	eststo sumstats10_`SUBSET':estpost tabstat hh_N_wtaxpayers 		[aw=_psweight] if $ifstatement		, statistics(N mean p25 p50 p75 )  columns(statistics)
}
restore

*[TABLE A.10] Propensity Score Weighting Robustness
foreach i of numlist 1(1)10 {
	#delimit;
	if `i' == 1  {; local decimals = 2; };
	if `i' == 2  {; local decimals = 0; };
	if `i' == 3  {; local decimals = 2; };
	if `i' == 4  {; local decimals = 0; };
	if `i' == 5  {; local decimals = 2; };
	if `i' == 6  {; local decimals = 2; };
	if `i' == 7  {; local decimals = 2; };
	if `i' == 8  {; local decimals = 0; };
	if `i' == 9  {; local decimals = 2; };
	if `i' == 10 {; local decimals = 2; };
	#delimit cr
	* count(fmt(%12.0fc)) dropped
	local cells = "mean(fmt(%12.`decimals'fc)) p25(fmt(%12.`decimals'fc)) p50(fmt(%12.`decimals'fc)) p75(fmt(%12.`decimals'fc))"
	#delimit; 
	esttab  sumstats`i'_0 sumstats`i'_1,			
			cells("`cells'") 
			substitute(	hh_gaver_tot_dummy 	"1[Giving$ _{h,2012} >0$]"
					hh_gaver_tot_w		"Giving$ _{h,2012}$ if $ >0$"
					hh_wtax_dummy   	"$ 1[wtax_{h,2012}>0]$"
					
					hh_wtax_b12_condMVHS    "$ wtax_{h,t}$ if MVHS$ _{h,2012}>0$"
					hh_wtax_b12		"$ wtax_{h,t}$ if $ >0$"
					hh_MVH_priseccabM_2012 	"MVH$ _{h,2012}$, MNOK"
					hh_MVH_secM_2012	"MVHS$ _{h,2012}$, MNOK"
					age			"Age$ _{i,2012}$"
					hh_MNW_MNOK			"Net Wealth$ _{h,2012}$, MNOK"
					hh_grossinc_MNOK		"Gross income$ _{h,2012}$, MNOK"
					hh_MWTR			"MWTR$ _{h,2012}$"
					hh_AWTR			"AWTR$ _{h,2012}$"
					hh_N_wtaxpayers		"Number of adults$ _h$"
					)	
			booktabs nolines nonum nomtitle collabels(none) noobs frag 
			posthead()
			prefoot(); 
	#delimit cr
}

* Run event regression with propensity score weighting
	disp "......... $S_TIME"
	#delimit;
	reghdfe	hh_logk3371
			year#c.hh_MVH_secM_2012#c.not2012
			c.hh_MVH_secM_2012 
			
		if $if_TVHS_cTNW_2 & hh_head_2012==1
		& inrange(year,${year_start},${year_end})
		[pw=_psweight]
		, cluster(hhid_2012)
		absorb(year##c.(	$base_control	
					hh_TNW_2012	hh_TNW_2012_sq 	hh_TNW_2012_cu
					age_2012 age2_2012 age3_2012 famsize_2012 famsize2_2012
					${inc_control}
					hh_N_wtaxpayers_2012
					hh_MVH_cabsec_dummy_2012
					
				)
			
			)
		;
	#delimit cr
	disp "......... $S_TIME"
	estimates save ${temp_estimates}/robust_event_PropensityScoreWeighting, replace
	mat TPSW = r(table)'

	
sum _ps if year==2012 & $if_TVHS_cTNW_2 & hh_head_2012==1 & hh_MVH_secM_d,d
	
*=================================================================================*

*=================================================================================*
capture drop _gps*
gen double _gps_cutpoints = .
replace _gps_cutpoints = 0 if hh_MVH_secM_2012 == 0
replace _gps_cutpoints = 1 if hh_MVH_secM_2012 > 0.25 	& !missing(hh_MVH_secM_2012)
replace _gps_cutpoints = 2 if hh_MVH_secM_2012 > 0.50 	& !missing(hh_MVH_secM_2012)
replace _gps_cutpoints = 3 if hh_MVH_secM_2012 > 0.75	& !missing(hh_MVH_secM_2012)
replace _gps_cutpoints = 4 if hh_MVH_secM_2012 > 1.00	& !missing(hh_MVH_secM_2012)
replace _gps_cutpoints = 5 if hh_MVH_secM_2012 > 1.50	& !missing(hh_MVH_secM_2012)
replace _gps_cutpoints = 6 if hh_MVH_secM_2012 > 2.00	& !missing(hh_MVH_secM_2012)
replace _gps_cutpoints = 7 if hh_MVH_secM_2012 > 2.50	& !missing(hh_MVH_secM_2012)
replace _gps_cutpoints = 8 if hh_MVH_secM_2012 > 3.50	& !missing(hh_MVH_secM_2012)
replace _gps_cutpoints = 9 if hh_MVH_secM_2012 > 5.00	& !missing(hh_MVH_secM_2012)
replace _gps_cutpoints = 10 if hh_MVH_secM_2012 > 7.00	& !missing(hh_MVH_secM_2012)

disp ".... $S_TIME "
#delimit;
capture gpscore	
			hh_MVH_priseccabM_2012
			hh_TNW_2012	hh_TNW_2012_sq 	hh_TNW_2012_cu
			age_2012 age2_2012 age3_2012 famsize_2012 famsize2_2012
			hh_loggrossinc_2012
			hh_N_wtaxpayers_2012
			if year==2012 & $if_TVHS_cTNW_2 & hh_head_2012==1
			,
		t(hh_MVH_secM_2012)
		gpscore(_gps_gpscore)
		predict(_gps_predict)
		sigma(_gps_sigma)
		cutpoints(_gps_cutpoints)
		index("mean")
		nq_gps(11);
#delimit cr
disp ".... $S_TIME "

capture drop _gps
by LOPNR (year): egen _gps = total(_gps_gpscore), mi

* Regression weights
/* see Hirano and Imens (2001): Estimation of Causal Effects Using [...] */
capture drop _gpsweight
gen 	_gpsweight = 1/_gps 		if hh_MVH_sec_dummy_2012==1 & !missing(_gps)
replace _gpsweight = 1/(1-_gps) 	if hh_MVH_sec_dummy_2012==0 & !missing(_gps)
sum _gpsweight,d
replace _gpsweight = . if _gpsweight > `r(p95)'
* Run event regression with GENERALIZED propensity score weighting
	disp "......... $S_TIME"
	#delimit;
	reghdfe	hh_logk3371
			year#c.hh_MVH_secM_2012#c.not2012
			c.hh_MVH_secM_2012 
			
		if $if_TVHS_cTNW_2 & hh_head_2012==1
		& inrange(year,${year_start},${year_end})
		[pw=_gpsweight]
		, cluster(hhid_2012)
		absorb(year##c.(	$base_control	
					hh_TNW_2012	hh_TNW_2012_sq 	hh_TNW_2012_cu
					age_2012 age2_2012 age3_2012 famsize_2012 famsize2_2012
					${inc_control}
					hh_N_wtaxpayers_2012
					hh_MVH_cabsec_dummy_2012
					
				)
			
			)
		;
	#delimit cr
	disp "......... $S_TIME"
	estimates save ${temp_estimates}/robust_event_GPropensityScoreWeighting, replace
	mat TGPSW = r(table)'

	
* Event plots

mat SPSW  = TPSW[1..13,1..6]
mat SGPSW = TGPSW[1..13,1..6]
estimates use ${temp_estimates}/event_hh_logk3371
estimates replay
mat Thh_logk3371 = r(table)'
mat Shh_logk3371  = Thh_logk3371[1..13,1..6]

foreach X in "hh_logk3371" "PSW" "GPSW"  {
	capture drop S`X'*
	svmat2  S`X', rnames(S`X'_rn) 
		
		replace  S`X'_rn =   substr(S`X'_rn,1,4)
		replace  S`X'_rn = subinstr(S`X'_rn,".","",.)
		replace  S`X'_rn = subinstr(S`X'_rn,"_","",.)
		replace  S`X'_rn = subinstr(S`X'_rn,"b","",.)
		replace  S`X'_rn = subinstr(S`X'_rn,"o","",.)
		gen long S`X'_xb =     real(S`X'_rn)
		drop S`X'2 S`X'3  S`X'4 
		
		gen S`X'_xb_shifted = S`X'_xb - 0.15
		
		replace S`X'5 = S`X'1 if S`X'_xb==2012
		replace S`X'6 = S`X'1 if S`X'_xb==2012
}
replace SGPSW_xb_shifted = SGPSW_xb + 0.15

*[FIGURE A.15] Robustness: Propensity Score Weighting
#delimit;
twoway  
	(connected 	    Shh_logk33711 	Shh_logk3371_xb if Shh_logk3371_xb >=2005.5,		lcolor(orange)    yaxis(1) msymbol(o) mcolor(orange%75) lwidth(medthick))
	(rcap Shh_logk33716 Shh_logk33715 	Shh_logk3371_xb if Shh_logk3371_xb >=2005.5,		color(orange%40) yaxis(1) lpattern(solid))
	
	(connected 	SPSW1 	SPSW_xb_shifted 	if SPSW_xb >=2005.5,		lcolor(blue)    yaxis(1) msymbol(s) mcolor(blue%75) lwidth(medthick))
	(rcap 	  SPSW6 SPSW5 	SPSW_xb_shifted 	if SPSW_xb >=2005.5,		color(blue%25) yaxis(1) lpattern(solid))
	(connected 	SGPSW1 	SGPSW_xb_shifted 	if SGPSW_xb >=2005.5,		lcolor(green)    yaxis(1) msymbol(d) mcolor(green%75) lwidth(medthick))
	(rcap 	 SGPSW6 SGPSW5 	SGPSW_xb_shifted 	if SGPSW_xb >=2005.5,		color(green%25) yaxis(1) lpattern(solid))
	
	, graphregion(color(white))
	ylabel(-0.01(0.005)0.005, axis(1) gmin gmax glcolor(gs14))
	
	xlabel(2007(2)2018)
	xmlabel(2007(1)2018, nolabel)
	ytitle("log of Charitable Giving", axis(1) color(black))
	xline(2012, lcolor(gs4%15) lwidth(vvvthick))
	ytitle("log of Charitable Giving (2012=0)", axis(1) color(black))
	legend(order(1 "Main specification" 3 "P-score weighting"  5 "Generalized p-score weighting") ring(0) pos(2) region(lcolor(gs12)) margin(r=8) bmargin(1 0.5 0.5 4) cols(1))
	
	xtitle("Year")
	xsize(1.75) ysize(1);
#delimit cr

graph save ${figures}/wdon1_robust_event_pscore_gpscore.gph, replace
graph export ${figures}/wdon1_robust_event_pscore_gpscore.pdf, as(pdf) replace


*=================================================================================*
* Robustness wrt controls, polynomials, etc 
*=================================================================================*

global Yvar = "hh_logk3371"

gen hh_TNW_2012_4 = (hh_TNW_2012/1000000)^4
gen age4_2012  = (age_2012/50)^4
gen famsize3_2012 = famsize_2012^3
gen hh_MVH_priseccabM_2012_sq = (hh_MVH_priseccabM_2012/1000000)^2
capture drop one
gen byte one = 1


global controls0 	= "one"
global controls1 	= "hh_MVH_priseccabM_2012"
global controls2 	= "hh_MVH_priseccabM_2012 			hh_TNW_2012	"
global controls3 	= "hh_MVH_priseccabM_2012 			hh_TNW_2012	hh_TNW_2012_sq 	hh_TNW_2012_cu	"
global controls4 	= "hh_MVH_priseccabM_2012 			hh_TNW_2012	hh_TNW_2012_sq 	hh_TNW_2012_cu			age_2012 "
global controls5 	= "hh_MVH_priseccabM_2012 			hh_TNW_2012	hh_TNW_2012_sq 	hh_TNW_2012_cu			age_2012 age2_2012 age3_2012"
global controls6 	= "hh_MVH_priseccabM_2012 			hh_TNW_2012	hh_TNW_2012_sq 	hh_TNW_2012_cu			age_2012 age2_2012 age3_2012 		 famsize_2012"
global controls7 	= "hh_MVH_priseccabM_2012 			hh_TNW_2012	hh_TNW_2012_sq 	hh_TNW_2012_cu			age_2012 age2_2012 age3_2012 		 famsize_2012  famsize2_2012"
global controls8 	= "hh_MVH_priseccabM_2012 			hh_TNW_2012	hh_TNW_2012_sq 	hh_TNW_2012_cu			age_2012 age2_2012 age3_2012 		 famsize_2012  famsize2_2012 			hh_loggrossinc_2012"
global controls9 	= "hh_MVH_priseccabM_2012 			hh_TNW_2012	hh_TNW_2012_sq 	hh_TNW_2012_cu			age_2012 age2_2012 age3_2012 		 famsize_2012  famsize2_2012 			hh_loggrossinc_2012 hh_N_wtaxpayers_2012"
global controls10	= "hh_MVH_priseccabM_2012 			hh_TNW_2012	hh_TNW_2012_sq 	hh_TNW_2012_cu			age_2012 age2_2012 age3_2012 		 famsize_2012  famsize2_2012 			hh_loggrossinc_2012 hh_N_wtaxpayers_2012 hh_MVH_cabsec_dummy_2012"
global controls11 	= "hh_MVH_priseccabM_2012 			hh_TNW_2012	hh_TNW_2012_sq 	hh_TNW_2012_cu	hh_TNW_2012_4	age_2012 age2_2012 age3_2012 		 famsize_2012  famsize2_2012 			hh_loggrossinc_2012 hh_N_wtaxpayers_2012 hh_MVH_cabsec_dummy_2012"
global controls12 	= "hh_MVH_priseccabM_2012 			hh_TNW_2012	hh_TNW_2012_sq 	hh_TNW_2012_cu	hh_TNW_2012_4	age_2012 age2_2012 age3_2012 age4_2012	 famsize_2012  famsize2_2012 			hh_loggrossinc_2012 hh_N_wtaxpayers_2012 hh_MVH_cabsec_dummy_2012"
global controls13 	= "hh_MVH_priseccabM_2012 			hh_TNW_2012	hh_TNW_2012_sq 	hh_TNW_2012_cu	hh_TNW_2012_4	age_2012 age2_2012 age3_2012 age4_2012	 famsize_2012  famsize2_2012 famsize3_2012	hh_loggrossinc_2012 hh_N_wtaxpayers_2012 hh_MVH_cabsec_dummy_2012"
global controls14 	= "hh_MVH_priseccabM_2012 hh_MVH_priseccabM_2012_sq	hh_TNW_2012	hh_TNW_2012_sq 	hh_TNW_2012_cu	hh_TNW_2012_4	age_2012 age2_2012 age3_2012 age4_2012	 famsize_2012  famsize2_2012 famsize3_2012	hh_loggrossinc_2012 hh_N_wtaxpayers_2012 hh_MVH_cabsec_dummy_2012"
global controlsALL 	= "hh_MVH_priseccabM_2012 			hh_TNW_2012	hh_TNW_2012_sq 	hh_TNW_2012_cu 			age_2012 age2_2012 age3_2012		 famsize_2012  famsize2_2012 			hh_loggrossinc_2012 hh_N_wtaxpayers_2012 hh_MVH_cabsec_dummy_2012"
	
foreach c of numlist 0(1)14{
	disp "================================================================"
	disp "======= $S_TIME:	c = `c' 			  ============"
	disp "================================================================"
	local controls = "${controls`c'}"
	disp "`controls'"
	
	disp "......... $S_TIME"
	#delimit;
	eststo stepwise_`c':
	reghdfe	$Yvar
			year#c.hh_MVH_secM_2012#c.not2012
				c.hh_MVH_secM_2012 
				
		if $if_TVHS_cTNW_2 & hh_head_2012==1
		& inlist(year,2012,2018)
		, cluster(hhid_2012)
		absorb(year##c.(
				`controls'	
			)
		)
		;
	#delimit cr
	estimates save ${temp_estimates}/robust_stepwise_controls_`c', replace
	
}

*[TABLE A.11]
#delimit; 
esttab  stepwise_*, 
	b(3) se(3)
	varwidth(30) stats(r2 N, fmt(2 0 )) 
	/*star(* 0.05 ** 0.01 *** 0.001)*/ nostar  nocons
	substitute (  	hh_MVH_secM_2012 "$ MVHS_{h,2012}$"
			2018.year#c.	"[2018]"
			#c.not2012	"$ 1[t\neq 2012]"
			R2		"R$ ^2$"
					
	)
	booktabs nolines nonum nomtitle  frag replace;
#delimit cr
	
